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Abstract 


Pulsar wind nebulae are fascinating systems, and archetypal sources for high-energy astrophysics in general. Due to 
their vicinity, brightness, to the fact that they shine at multi-wavelengths, and especially to their long-living emission 
at gamma-rays, modelling their properties is particularly important for the correct interpretation of the visible 
Galaxy. A complication in this respect is the variety of properties and morphologies they show at different ages. 
Here we discuss the differences among the evolutionary phases of pulsar wind nebulae, how they have been modeled 
in the past and what progresses have been recently made. We approach the discussion from a phenomenological, 
theoretical (especially numerical) and observational point of view, with particular attention to the most recent results 


and open questions about the physics of such intriguing sources. 
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1 INTRODUCTION 


The violent death of a massive star (M z 8Mọ) as a supernova 
(SN) is thought to leave behind a compact remnant, in many 
cases in the form of a rotating and magnetized neutron star, 
known as pulsar (PSR). Pulsars have outflows, blowing out 
from the open region of their magnetosphere in the form of 
highly relativistic (with Lorentz factor I, >> 1), magnetized 
(with ratio between the Poynting and kinetic energy fluxes 
o > 1) and cold plasma winds. This wind blows inside 
the cold debris of the parent star (the ejecta), themselves 
slowly expanding (with typical speeds ~ 300 — 5000 km 
s^, Chevalier 1976) in the outer inter-stellar medium (ISM). 
The interaction between the relativistic pulsar wind and these 
confining ejecta, gives rise to the formation of a wind bubble 
bounded in the inside by a strong magneto-hydrodynamic 
(MHD) termination shock (TS). This literally terminates the 
pulsar wind: at its surface the plasma is heated due to the 
randomization of the bulk motion, magnetic field is dissipated, 
and particles are accelerated (Gaensler & Slane, 2006; Slane, 
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2017). This scenario typically applies to young pulsar wind 
nebulae but, as we will see in the following, a rather similar 
one can explain the formation of bow shock nebulae created 
by the outflow emanating from an (older) pulsar directly 
interacting with the ISM. 


The pulsar wind nebula (PWN) becomes observable due 
to the non-thermal emission of particles accelerated at the 
shock. For young systems most of the energy is radiated as 
synchrotron emission from the relativistic e^ — e* pairs spi- 
ralling in the local magnetic field. Higher energy gamma-rays, 
in the GeV to TeV ranges, are instead produced by inverse 
Compton scattering (IC) processes between the same relativis- 
tic particles responsible for synchrotron, and various back- 
ground photons (either synchrotron photons, photons from 
the infrared background, photons of the cosmic microwave 
background, starlight photons). The high-energy spectrum 
might extend up to hundreds of TeV, where Klein-Nishina 
suppression produces a natural cut-off. The appearance of old 
systems, as we will detail in the following, might be very dif- 
ferent, both due to lower injection of energy from the pulsar, 
the synchrotron cooling and the value of the nebular mag- 
netic field (with IC possibly even exceeding the synchrotron 
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Figure 1. Sketch of the structure of a PWN: the nebula (in light red) is 
embedded in the expanding SNR ejecta (in aquamarine color), separated 
from the outer ISM by the SN forward shock. The darker region at the PWN 
center mimics the pulsar wind termination shock, a region producing no 
emission and thus visible as an under-luminous area in the inner nebula. The 
inner structure of the PWN, characterized by a jet-torus morphology visible 
at X-rays, is drawn in light blue. 


emission). 

Despite the success of purely leptonic models in reproduc- 
ing the spectral energy distribution of many PWNe, room 
is left for the possible presence of hadronic emission, pro- 
duced in the decay of neutral pions following proton-proton 
collisions. The presence of hadrons in the pulsar wind is still 
matter of discussion (Atoyan & Aharonian, 1996; Bednarek 
& Protheroe, 1997; Amato et al., 2003; Bednarek & Bartosik, 
2003), and a clear evidence for their existence, or proof of 
their absence, is missing. For example, the recent detection 
by the LHAASO experiment of photons with energy above 
1 PeV (i.e. 10? eV) in coincidence with the position of the 
Crab nebula (Cao et al., 2021a), has been claimed as a pos- 
sible indication of high energy protons in the Crab wind. 
Unfortunately, despite this very impressive result, the lack of 
a robust statistics above 500 TeV makes the overall spectrum 
still compatible with a fully leptonic scenario (for a more in 
depth discussion see Amato & Olmi 2021). 

Much of the importance of PWNe is tied to the fact that 
the Crab nebula, the prototype of this class, is the brightest 
non-thermal object in the sky over almost its entire electro- 
magnetic spectrum, from radio, to optical, X-rays and beyond. 
As such, it offers us a unique laboratory where high energy 
astrophysical processes can be investigated in great details. 
On the other hand, the fact that this class includes many other 
systems, even if fainter and not as well constrained, makes it 
possible to generalize many of the findings. From the morpho- 
logical point of view, PWNe appear as fill-center objects, with 
possibly a darker region surrounding the pulsar in the very 
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inner nebula that marks the pulsar wind zone (the dark blue 
region in the sketch of Fig. 1). Due to its physical properties, 
the pulsar wind does not produce any sizeable amount of radi- 
ation (apart perhaps pulsed emission), and it is then extremely 
hard to infer information on the pulsar outflow properties be- 
fore it crosses the TS. In the Crab nebula, the extension of 
the emitting region decreases with increasing frequency from 
radio to X-rays as a direct effect of the synchrotron process: 
the particles life-time against synchrotron losses is in fact 
inversely proportional to their energy. Radio emission then 
appears in general more extended (and uniform) than the 
X-ray one, that is mostly confined to the inner nebula, and 
highlights the complex structure of the local plasma. A very 
rough comparison of the difference in extension can be seen 
in the red sketch of Fig. 1 (light red for the radio nebula, light 
blue for the X-ray one). This however is not the case in all 
the other systems. In MSH15-52 the radio, X-ray and gamma- 
ray sizes are comparable (Gaensler et al., 2002; Aharonian 
et al., 2005; Leung & Ng, 2016). In G21.5-0.9 the infrared 
(IR) size is smaller than the X-ray one, that indeed is very 
similar to the radio extension (Guest et al., 2020; Zajczyk 
et al., 2012). In the Kookaburra PWN the X-ray size is larger 
than the radio one (Roberts et al., 1999, 2001). This just to 
indicate that extrapolating global trends from the Crab nebula 
might be misleading and that there is a large diversity in the 
population. 

To date around 110! sources, observed in different energy 
bands, have been recognized as — or candidate to be (~ 20 
cases) - PWNe. Around 20 of those systems do not have an 
associated pulsar. Most of them have been detected first at 
X-rays, while the multi-wavelength identification of a PWN 
is, in general, not an easy task. The extended radio emis- 
sion might be contaminated by diffuse emission from the 
surroundings, and from radio data it is very difficult to con- 
strain precisely the morphology of the system, which can help 
to identify it. A more important reason of this difficulty is the 
lack of sensitivity to the required large angular scales (typical 
of old PWNe) of most radio interferometers which do these 
observations. Optimal strategies would require the combined 
use of interferometry and single-dish observations (Kurono 
et al., 2009). X-ray emission is then typically the main dis- 
criminant for the identification, especially for young PWNe, 
given that there are not so many other extended Galactic 
sources bright in this band. Since the beginning of the cen- 
tury, mainly thanks to the very impressive sensitivity of the 
Chandra telescope for X-ray imaging, our ability to derive 
information about PWNe morphology has improved signifi- 
cantly. Unfortunately, X-rays are the first to fade away as time 
passes by: X-ray emission only lasts for a (relatively) small 
fraction of the PWN life. Middle aged PWNe show in fact 
very limited (if not any) X-ray emission, making them hardly 
detectable. On top of this, as the PWN gets larger, issues with 
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the limited field of view of many of the X-ray instruments 
begin to play a role. 

When getting older, PWNe end up being observable mostly 
at gamma-rays, where the dominant population of emitting 
particles is the same producing the long lasting radio emis- 
sion in the synchrotron band. Unfortunately, the instruments 
resolution at gamma-rays is much worse than that at lower 
energies, and it is then very difficult to infer the nature of a 
source from its morphology. Present gamma-ray data almost 
certainly contain a larger number of unidentified PWNe than 
of identified ones: 14 are the PWNe firmly identified in the 
last H. E. S. S. Galactic Plane Survey (Abdalla et al., 20182), 
while around 45 are the unidentified sources. Many of these 
— if not almost all — are actually believed to be unidentified 
PWNe, expected to represent the most numerous class of 
sources in the very high energy sky (de Ofia-Wilhelmi et al., 
2013; Klepser et al., 2013; Abdalla et al., 2018b). 

In the next future we will see the advent of a new gen- 
eration of Imaging Atmospheric Telescopes (IACTs) for 
gamma-ray observations, as the Cherenkov Telescope Array 
(CTA, Actis et al. 2011) or the ASTRP Mini-Array (Scuderi 
et al., 2022). If compared with water Cherenkov detectors, as 
LHAASO or Tibet As-y, they will operate in a reduced energy 
range but with a much better sensitivity and resolution. As a 
consequence, the number of PWNe detected at gamma-rays 
is expected to increase significantly (Fiori et al., 2022; Remy 
et al., 2022), posing the problem of the identification of a 
very large number of new PWNe, mostly in their middle-age 
stage. 

In this review we will discuss the different evolutionary 
phases of a PWN, from its early quasi-spherical evolution, to 
the late phases characterized by a completely different shape 
and possibly by strong asymmetries, presenting available the- 
oretical models, results form numerical studies, observational 
hints and open problems. 

The review is structured as follows: in Section 2 we qual- 
itatively introduce the properties of the different phases in 
which we can roughly divide the evolution of a PWN. Here 
we also consider PWNe evolving in different environments, 
discussing their nature as prototype high-energy astrophysi- 
cal sources. In the following Section 3, we will briefly inspect 
theoretical and numerical models of PWNe from an historical 
point of view; we will highlight results from recent 3D MHD 
models of young and old PWNe. In Section 4 we describe 
the radiation processes at the base of the observed multi- 
wavelength emission from those sources, and discuss the 
possible acceleration mechanisms producing particles respon- 
sible for the observed emission. A description of the actual 
sample of PWNe, from a multi-wavelength point of view,is 
presented in Section 5, where we also review the observa- 
tional properties of systems at different stages of evolution. 
Finally, in Section 6 we discuss the future prospects both in 
terms of observations and modelling, with particular focus 
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on the actions needed to answer old and new open questions 
in pulsar and PWNe physics. 

A brief summary of the main arguments treated in this review 
and our concluding remarks can be found in Section 7. 


2 EVOLUTIONARY STAGES OF PWNE 


The evolution of a PWN can be ideally divided into four 
different stages, even if a clear distinction to mark the onset of 
a new phase from the end of the previous one is hard to make 
especially for the later stages, which approach a continuum. 
A graphical representation of the evolutionary path of a PWN 
— and the four phases — can be seen in Fig. 2. The typical 
duration of each phase depends on both the evolution of the 
pulsar properties (mainly the spin-down luminosity) and the 
dynamics of the supernova remnant (SNR) where the PWN 
expands. 


2.1 PSR Evolution 


The engine driving the dynamics of the PWN is the central 
pulsar. It loses energy, injecting it in the form of a relativistic 
plasma, following a typical magnetic dipole spin-down rate 
(or luminosity, see Pacini & Salvati 1973): 


: P 
E-L-4mIa, (1) 


where J ~ 109 g cm? is the pulsar moment of inertia, while P 
is its rotational period and P its time derivative. For canonical 
pulsars the period varies typically between 0.01 s and 1 s, 
while the period derivative is in the range 107? — 107!! ss~}, 
As one might expect, X-ray bright synchrotron nebulae, likely 
younger, tend to be associated with energetic PSRs (E = 
1096 erg s~!). Older objects instead tend to be associated with 
less powerful engines (see e.g. Fig. 8 in the following Sec. 5). 

To model the evolution of the pulsar energy injection with 
time, it is customary to assume that its spin slows down from 
an initial value Po according to: P œ P?-", where the coeffi- 
cient n is called the braking index. This is usually assumed to 
be in the range 2 < n € 3, with n = 3 the case of pure dipole 
spin-down. However, recently larger values have also been 
suggested (as in the case of H. E. S. S. J1640-465, reported 
in Archibald et al. 2016), even if their physical meaning is 
still not fully understood (Parthasarathy et al., 2020). If n 
remains constant during the pulsar life, the variation of the 
pulsar energy input is (Pacini & Salvati, 1973): 


t -(n*1)/(n-1) 
, (2) 


L(t) = Lol + — 
To 


where To = Po/[(n— 1)Po] is the initial spin-down time of the 
pulsar. From this equation it is easy to see that the pulsar input 
can be considered as constant only for t « To, while properly 
accounting for its temporal evolution becomes important for 
the correct modelling of the nebula beyond Tọ. 
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OUT of REVERBERATION OUT of SNR 


Figure 2. Sketch of the different evolutionary phases of a PWN. From left to right: the PWN expands in the un-shocked ejecta; the RS crushes the PWN and 
marks the begin of reverberation; the outcome of reverberation depends on the possibility for the PWN to efficiently contrast the compression exerted by 
SNR; high-velocity enough pulsars eventually exit their parent SNR bubble and interact directly with the ambient medium, producing cometary nebulae (bow 


shocks). The colors for the PWN ant the ejecta are the same used in Fig. 1. 


22, SNR Evolution 


A very important point that we want to make before dis- 
cussing the phenomenology of the various evolutionary 
phases, is that the age of the PWN (or the pulsar) is not 
a good indicator of the evolutionary stage of the PWN. A 
much better indication is instead given by the SNR charac- 
teristic time, first introduced by Truelove & McKee (1999), 
providing the typical time scale for the evolution of the SNR: 
= 5/6 -1/3 

fch = E PMY oe 2 (3) 
where E,, is the SN explosion energy (typically of order 
of 10°! erg), Mej is the mass of the ejecta (in the range ~ 
6 — 20M, for SNRs hosting a PWN, see e.g. Smartt 2009) 
and pism the density of the interstellar medium (ISM). 

We wish to remark here that, contrary to the widely used 
Truelove & McKee (1999) model, which predicts large varia- 
tions, depending on the progenitor structure, in the evolution 
of the SNR reverse shock (RS, an important factor in set- 
ting the various PWN evolutionary stages), and in particular 
on the time it takes to reach the center of the remnant (an 
event that we name implosion, happening at time fimpio), a 
much smoother trend has been recently found by Bandiera 
et al. (2021), using fully Lagrangian simulations. In particular, 
large variations of the implosion time only appear when dras- 
tically chancing the structure of the ejecta core (see e.g. Fig. 4 
of that work), while for the largely used case in PWN+SNR 
models (ejecta with a flat core plus a steep envelope), implo- 
sion happens at timplo = 2.4 ten. This value clearly represents 
the maximum time the PWN can remain in free-expansion 
(see next subsection), and can be used as a sort of theoretical 
upper limit for the duration of this first phase see Sec. 2.4. 
In reality, for typical PSR energy injection, the duration of 


the free-expansion phase is no longer than about half of the 
implosion time. The implosion time is more representative 
of the time at which the compression due to the interaction 
of the PWN with the SNR reaches the maximum. It should 
be clear that, depending on the properties of the SNR, this 
time might be very different for diverse systems, and this is 
the reason why the age of the PWN (or the pulsar) is not a 
good indicator of the stage of its evolution. Indeed, a change 
by a factor of 2 in the mass of the ejecta leads to a 1.8 factor 
of difference in the implosion time: if we assume, for exam- 
ple, a fixed SN energy of Es = 10°! erg, and an ISM with 
number density 1 particle cm~?, a PWN in a 6Mo SNR could 
stay at most for ~ 5000 yr in the first stage, while a PWN 
powered by the same pulsar in a 10M, remnant could stay in 
free-expansion for ~ 9000 yr. 

Bandiera et al. (2023) show that the entire evolution of 
each PWN-SNR system is fully determined — in the ab- 
sence of significative radiative losses — by the two quantities: 
[To/tch, Lotch/Esn]. These in practice weight the pulsar time 
and energetics with respect to the SNR ones. 

In the following of this section we limit ourselves to a 
phenomenological discussion of the properties of the various 
phases. 


2.3 Free-expansion 


In the first phase the PWN expands, with a mild acceleration, 
in the cold freely-expanding ejecta of the SNR core, hence 
its name. The typical PWN expansion speed is of the order of 
few thousands km s~!, much higher than the average velocity 
the PSR can acquire during the SN event (the kick velocity, 
ranging in 100-500 km s^, see e.g. Faucher-Giguére & Kaspi 
2006), so that one can safely consider the PSR to be station- 
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ary at this stage. During this phase the evolution of the PWN 
is independent from that of the SNR shell, since no interac- 
tion has been established yet between the two (Reynolds & 
Chevalier, 1984). Based on the results by Jun (1998), showing 
that the PWN collects material in a thin-shell at its boundary 
during its initial expansion (t « To), a simplified description 
of the PWN evolution can be obtained using the following 
formulas: 


d 
7 UnPORO') = LORO, (4) 
d dR(t)\ | 2, dM() RO 

zl Dd = 4x POR +- y 0) 


known as the “thin-shell approximation". Here P(t) is the 
PWN pressure, M(t) the mass of the thin-shell and R(t) the 
shell radius, that within the thin-shell approximation, can be 
taken as that of the PWN. 

The freely-expanding ejecta have a density profile char- 
acterized by a flat core surrounded by a steep envelope. It 
is customary to use power laws to describe them: the steep 
envelope as r^^? (with w > 5) and the shallow core as r^? 
(with 6 « 3, see Bandiera et al. 2021 and references therein): 


A (u/rf 0, ifr « ut, 
PCO ou : (6) 
A(u/r)'t^^, ifr ut, 


with the parameters A and v (the expansion velocity of the 
ejecta core) that depend on the SN energy and mass of the 
ejecta as: 


(5 a ôw a 5) Esn 
2m(w-6) yp’ M 


|. RG- 8Xe - 5) E; 
* = A G-3-—3) Ma ` 


(8) 


For the simplified case of the flat density profile (6 = 0) 
plus a steep envelope (w = oo), actually the most commonly 
used, the mass of the shell can be expressed as: M(t) = 
4nR?(t)A/(3P). At early enough times, when the pulsar input 
can be still considered as constant (namely t « To, so that 
L(t) = Lo) an analytical solution for Eq. 4-5 can be found 
(Bandiera et al., 2023): 


E (3-5 6) Lo 16-2 1/(5-6) 
rem | (9 — 25)(11 — 26) 4n Ai? 


Rpwn(t) (9) 


This reduces to the well known trend found for the specific 
case of a flat core (6 = 0) by Reynolds & Chevalier (1984) 
and later by van der Swaluw et al. (2001): Rpwn(a) œ p 
On the contrary an analytic solution valid at all times up to 
reverberation has not been found, while a tentative solution 
based on the series expansion of the spin-down law was 
investigated in Bucciantini et al. (2003). For the flat core 
case, Bandiera et al. (2023) have recently shown that a simple 
solution can be determined based on the fitting of numerical 
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Figure 3. Left panel: time evolution of the PWN radius (in units of Voro 
and To) for four values of the braking index in the range n = (1.8; 4.1}. In 
magenta color we show the pure dipole braking (n = 3), while in orange the 
value measured for Crab (n = 2.33, Lyne et al. 2015; Horvath 2019). Right 
panel: ratio of the previous two. 


results: 
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Rows O|. = Voto (10) 


where the parameters c;, a and b depend on the value of the 
braking index n, while Voto = 1.91(LoTo | Eg) Poft) Ras, 
and Ron = M yon is the characteristic radius. A very in- 
teresting result is that the evolution of the PWN in the free- 
expansion phase is poorly affected by the value of the braking 
index, at least in the range 1.8 < n < 4.1. This can be seen in 
Fig. 3, comparing the radial evolution computed with Eq. 10 
for n = 1.8 (c; = 0.86310, a = 0.78492, b = 0.76490), 
n = 2.33 (c; = 0.9517, a = 0.73014, b = 0.71979), n = 3 
(c; = 1.0329, a = 0.66355, b = 0.65937) and n = 4.1 
(c; = 1.1334, a = 0.59129, b = 0.59210), up to a maxi- 
mum time of 5 — 1079, much longer than the duration of this 
initial phase for typical systems. 

As we will better see in Sec. 3, being characteristic 
of young objects, and hence of bright systems, the free- 
expansion phase has received a lot of attention in the past 
in terms of modelling, with a variety of approaches based 
on different strategies (Reynolds & Chevalier, 1984; van der 
Swaluw et al., 2001; Bucciantini et al., 2003; Gelfand et al., 
2009; Martín et al., 2012; Fiori et al., 2022; Komissarov & 
Lyubarsky, 2004; Del Zanna et al., 2006; Porth et al., 2014; 
Olmi et al., 2014, 2016). 


2.4 Reverberation 


When the PWN outer boundary hits the SNR RS, the free- 
expansion phase ends and reverberation starts. In reality, due 


to Rayleigh-Taylor like instabilities (R-T), that develop at the 
PWN boundary, as well as asymmetries of different nature 
(e.g. high proper motion of the pulsar, strong ISM density 
gradients, asymmetries induced from the stellar explosion), 
this transition is not as sharp as simplified 1D models would 
suggest, neither is the following reverberation phase. This 
phase and its transition have received much less attention due 
to several reasons: their extreme complexity, the fact that the 
SNR evolution and properties become more relevant, the PSR 
kick velocity begins to play a role. Moreover, there are few 
well characterized systems from an observational point of 
view that can be taken as benchmark to evaluate the accuracy 
of a model. Only recently a detailed study of the transition be- 
tween free-expansion and reverberation have been presented 
(Bandiera et al., 2023), still limited to a simplified 1D evolu- 
tion. The time at which this transition happens (that we name 
fbegrev) is technically obtained as the intercept between the RS 
radius and the PWN one. It depends on the relation between 
the PWN energetics (Loro) and Esn, and can be computed 
with some accuracy in the limit To < ten using the following 
formula by Bandiera et al. (2023), expressed for simplicity in 
terms of the variable Ag = log, )[(LoTo)/Esnl: 
— exp(—0.1494 + 1.1606 Ag) 


1 
ficte) = 2.4102 ta. (11 
egre (AE) 1 + exp(1.6831 + 0.680547) ^ PD 


An extension to larger ro is possible substituting to LoTo 
the energetic of a massive shell interacting with the SNR 
in place of the PWN, as discussed in Sec. 5.1 of Bandiera 
et al. (2023). The time fpegrey typically ranges between 1 and 
2 characteristic times for most of the systems (see Fig. 2 of 
Bandiera et al. 2023). 

When reverberation starts, the PWN begins to interact 
directly with the shocked ejecta in the SNR. The pressure ex- 
erted by the ejecta generally induces a strong deceleration of 
the shell accumulated at the PWN boundary (van der Swaluw 
et al., 2001) which, in many cases, turns into a compression 
of the nebula. During this compression the internal energy of 
the PWN, and its magnetic field, increase. Ultimately, depend- 
ing on the efficiency of radiation losses and magnetic energy 
dissipation, the total internal pressure rises to become com- 
parable to the external one, and the compression is suddenly 
reverted into a new expansion. 

One-zone and 1D models predict a number of successive 
oscillations and expansions before the system reaches a new 
equilibrium in the Sedov-Taylor phase, and the PWN keeps 
expanding without bouncing back anymore. This oscillatory 
behavior is what gives the name reverberation to this phase. 
As already pointed out in Bucciantini et al. (2011), this ex- 
tended oscillatory behavior is an artifact of the 1D approxi- 
mation. Higher dimensional models in fact show that during 
reverberation the PWN boundary is largely subjected to R-T 
like instabilities at the boundary (Blondin et al., 2001; Kolb 
et al., 2017), producing effective mixing between the PWN 
material and the outer one. This acts as a sort of viscous 
term that might halt oscillations already after the first com- 
pression. This is why in one-zone models reverberation is 
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usually stopped artificially after the first bounce (Martin et al., 
2012; Torres et al., 2014), or when the PWN pressure reaches 
equilibration with the outer one (Bucciantini et al., 2011). In 
simplified 1D models, the effect of the reverberation phase 
on the dynamics of the PWN can be measured in terms of 
the compression factor (CF), namely the ratio between its 
maximum radius (close to the begin of reverberation) and the 
radius at the maximum of compression, i.e the minimum one. 
More generally the parameter of interest is the change in the 
total volume taken by the relativistic non-thermal plasma of 
the PWN, since it translates in a variation of the nebular mag- 
netic field (B œ R^?). It is in fact this quantity that directly 
impacts the chance an old PWN might become visible again 
in X-rays. 

Roughly speaking, the effects of the reverberation can be 
divided into two extreme cases: (1) the PWN is sufficiently 
energetic, then the compression is relatively small, or even 
not appreciable, and the PWN continues its expansion almost 
undisturbed; (ii) the PWN is weak and then overwhelmed by 
the SNR pressure, contracting down to very small radii (vol- 
umes). Low energetic systems are the most critical in terms 
of modelling, since they might undergo violent compression 
with important modifications of their multi-wavelength spec- 
tral properties. In fact compression enhances the magnetic 
field and energises particles, leading to an increase in radiative 
losses. For very high compression efficiency (CF 2 1000), 
the PWN can enter a fast cooling regime, where a large frac- 
tion of the particle energy is lost, and the particle energy dis- 
tribution is strongly modified (Torres et al., 2019). Bandiera 
et al. (2023) show that these extreme behaviour is not ex- 
pected to be common within the present PWNe population, 
with only a relatively limited number of objects having pos- 
sible CF larger than few hundreds. On the contrary, the vast 
majority of the systems undergo a rather small compression, 
with CF ranging from a few to tens. However Bandiera et al. 
(2023) do not consider the effect of radiative losses and thus, 
this estimate might change depending on the relevance of 
losses in the first evolutionary phase. 


2.5 Post-Reverberation 


The reverberation phase, together with the pulsar proper mo- 
tion, is what shapes the PWN in the later evolutionary phases. 
Gradients in the ambient medium density are known to im- 
pact the evolution of the SNR reverse shock (Ferreira & de 
Jager, 2008; Kolb et al., 2017) and, in general, one expects 
that as time passes the level of asymmetries in the system 
grows. This implies that simplified one-zone or 1D models 
became progressively less accurate in the description of these 
systems, and attempts to include these extra effects (Gelfand 
et al., 2009) less and less reliable. On top of this, the mixing 
due to R-T instability can be so strong as to disrupt com- 
pletely the PWN as a coherent object. It is then clear that a 
PWN in this stage would likely be very far from being spheri- 
cally symmetric, rather having a complex distorted shape. As 
a consequence, modelling this phase can be very demanding: 
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the need to properly account for asymmetries requires the 
use of multi-D models. Being the shape and properties of 
the evolved PWN strongly dependent on its previous history, 
one has to follow its entire evolution, which implies a large 
dynamical range in term of temporal and spatial scales. More- 
over the morphology of the system, especially in the presence 
of mixing and instabilities, is very dependent on the model 
dimensionality and a comprehensive description can only be 
done in 3D. A very beautiful example of a 3D model of a 
largely asymmetric evolved system can be found in Kolb et al. 
(2017). 

Unfortunately there are not that many systems, with a de- 
tailed and robust characterization, that can be used to bench- 
mark our theoretical models for this evolutionary phase. Some 
of them like G327.1-1.1 (Temim et al., 2015; Ma et al., 2016), 
IC443 (Swartz et al., 2015), W44 (Frail et al., 1996; Petre 
et al., 2002), have been the targets of previous analyses (van 
der Swaluw et al., 2004; Bucciantini et al., 2011), but they 
are too few to provide a significative sample. 

From the one-zone description, one expects that when 
reverberation ends, with the oscillatory behaviour almost 
completely dissipated, the PWN relaxes to a steady subsonic 
expansion (van der Swaluw et al., 2001). This was found to 
happen on a long time-scale of ~ 20f,,. This time is possibly 
even longer, since the relaxation of the contact discontinuity 
of the SNR has been recently show to happens on times of 
~ 35ten, while between 20f,, and 30f, small oscillations are 
still present in the SNR radius (Bandiera et al., 2021). 


2.6 Late Phase - bow shocks 


Ultimately PWNe are likely to end their life as bow shock 
nebulae (BSPWN), due to the fact that a large fraction of 
pulsars is born with high kick velocity (100-500 km s“!, 
Cordes & Chernoff 1998; Arzoumanian et al. 2002; Faucher- 
Giguére & Kaspi 2006; Sartore et al. 2010; Verbunt et al. 
2017), and as such it is bound to emerge out of the parent 
SNR before the pulsar spin-down luminosity becomes so 
weak as to make particle acceleration and non-thermal syn- 
chrotron emission negligible. Considering the typical SNR 
decelerated expansion of the Sedov-Taylor phase, one can 
easily estimate in a few tens of thousand of years the time the 
pulsar takes to emerge the SNR bubble, to be compared with 
the longer typical age of pulsars in the Galaxy (^ 10° yr). 
From that moment on, the pulsar interacts directly with the 
ambient medium. Given the high speed of the pulsar, larger 
than the typical sound speed of the surrounding medium, its 
motion generally becomes supersonic already inside the SNR 
(Gaensler & Slane, 2006). IC443 and W44 are exemplary 
cases of BSPWNe still confined within their parent SNR 
(Swartz et al., 2015; Frail et al., 1996; van der Swaluw et al., 
2004). The supersonic motion induces the formation of a bow 
shock around the pulsar and its nebula, reshaping drastically 
the PWN. Now it does not appear as a bubble anymore but, 
on the contrary, it assumes an elongated cometary-like shape. 

In BSPWNe the pulsar is located at the bright head of 


a very elongated tail, extending in the direction opposite 
to the pulsar motion. The new morphology of the system 
depends on the balance of the ram pressure of the pulsar 
wind with respect to that of the incoming (in the frame of 
the PSR) ambient medium. This balance actually determined 
the thickness of the PWN at the head front, representing the 
characteristic dimension of a BSPWN, known as stand off 


distance: 
L 
dy = IL ÀÀ. (12) 
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where p is the density of the ambient medium. Gradients 
in the ambient medium reflect into asymmetries in the bow 
shock shape (Vigelius et al., 2007). On the contrary, even 
larger asymmetries in the pulsar wind energy distribution 
tend to have minor effects on the overall bow-shock shape 
(Olmi & Bucciantini, 20192), remaining mostly concentrated 
in the head, that is typically not resolved even in the deepest 
Chandra's observations. 

The PSR wind, shocked in the head, is now diverted along 
the tail, where the plasma becomes less magnetized and more 
turbulent with increasing distance from the pulsar (Wilkin, 
1996; Bucciantini & Bandiera, 2001; Bucciantini, 2002). How 
turbulent and magnetized the tail is depends not only on the 
PSR wind magnetization, but also on the level of asymmetry 
in its energy distribution, on the inclination of the PSR spin 
axis with respect to the direction of the kick velocity, and on 
the relative inclination and strength of the ambient medium 
magnetic field (Olmi & Bucciantini, 20192). 


2.7 Other environments 


PWNe represent proto-typical systems where one can wit- 
ness the interaction of a relativistic magnetized outflow with 
a confining surrounding environment. As such they form a 
paradigm for a variety of different classes of high-energy 
astrophysical objects. Moreover, PWNe can be found more 
or less anywhere a NS is active, even if the activity is just a 
transient one. Indeed, from an historical point of view, they 
have often been “the first ones" where high-energy astrophys- 
ical processes have been discovered/studied/understood, and 
much of the theory developed for their study has found appli- 
cation elsewhere. Here we briefly review some of these other 
environments. 


e PWNe, from the point of view of fluid dynamics, 
acceleration properties, and emission mechanisms, are 
representative of the very broad class of relativistic wind 
bubbles. The physics that have been thoroughly investi- 
gated to explain the observed presence of non-thermal 
particles or the properties of relativistic outflows at large 
distances from their engine (Lyubarskij, 1992; Li, 1993; 
Barkov & Komissarov, 2008; Chen & Zhang, 2021; 
Zhu et al., 2019; Kagan et al., 2015; Sironi et al., 2015) 
have proven to be quite general and with a much larger 
applicability, in terms of astrophysical systems ranging 


from gamma ray bursts (GRBs, Thompson, 1994) to 
active galactic nuclei jets and radio lobes (Uzdensky, 
2016, 2018). 


Bright persistent X-ray emitting PWNe require a steady 
engine, providing the necessary particle and energy 
injection. It was held for a long time that only canonical 
PSRs, with their active magnetospheres, capable of 
supporting pair creation cascades with high multiplicity 
(Timokhin & Harding, 2019) could be surrounded by 
such nebulae. However, recently has emerged a more 
dynamical picture, where transient PWNe can form and 
shine, even around neutron stars with inefficient pair 
creation, as in the case of magnetars and/or rotating 
radio transients (RRATs, Camero-Arranz et al., 2013; 
Vink & Bamba, 2009; Tanaka, 2016; Torres, 2017; 
Granot et al., 2017; Blumer et al., 2017; Torres, 2018; 
Margalit & Metzger, 2018). These transient nebulae, 
typically associated with bursting/active phases, are still 
electromagnetic powered, and shine in synchrotron, even 
if it is still unclear if the neutron star activity leads to 
freshly injected particles, or simply re-energizes particles 
already present and injected at earlier times. Even the 
driving energy reservoir is not well constrained in general, 
though magnetic energy has been suggested for the 
case of magnetars. Moreover, violent events, as giant 
flares, can produce transient nebulae, whose evolution 
resembles, even if at a faster pace, that of regular PWNe 
(Gaensler et al., 2005; Gelfand et al., 2005). 


PWNe can also form in binary systems, where the 
pulsar wind is confined by the ram pressure of the 
companion star. In case of massive stars, with strong 
equatorial flows, this interaction leads to the formation 
of a bow-shock, whose dynamics, and by consequence 
emission, can be both highly variable due to the possible 
large orbital eccentricity, and to the fact that it is also 
strongly affected by centrifugal and Coriolis forces. As in 
normal bow-shocks, particles can be accelerated to high 
energy, but unlike regular systems, the orbital dynamics 
can lead to a substantial mixing with the stellar material, 
a high level of turbulence, and the development of 
multiple shocks, with distributed acceleration. Moreover 
the presence of a bright source makes IC a relatively 
strong cooling process, to the point that it is unclear if the 
observed X-ray emission is synchrotron or IC (Neronov 
& Chernyakova, 2007; Kargaltsev et al., 2014; Bednarek 
& Sitarek, 2013; Zdziarski et al., 2010; Paredes-Fortuny 
et al., 2015; Bosch-Ramon et al., 2012). These systems 
offer us a unique opportunity to investigate the PSR 
wind and its acceleration properties, at distances much 
smaller than in regular PWNe, and can have important 
consequences on a large set of dissipative wind models 
(Kirk & Skjzraasen, 2003; Kirk, 2006; Lyubarsky, 2010; 
Takamoto et al., 2012; Cerutti & Philippov, 2017; Cerutti 
et al., 2020; Huber et al., 2021a,b). 
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e Proto-magnetar wind nebulae, and in general relativistic 
outflows from millisecond rotating newly born magnetars, 
have been invoked to explain both long (Usov, 1992; 
Thompson, 2007; Bucciantini et al., 2007, 2008; Metzger 
et al., 2011b) and short GRBs (Bucciantini et al., 2012; 
Wang, 2017; Metzger & Piro, 2014; Rezzolla & Kumar, 
2015), mostly because the detection of the so called 
“late activity" (Gompertz et al., 2014; Rowlinson et al., 
2013) has made almost mandatory to assume long lived 
engines, hardly compatible with the timescale for disk 
accretion in a stellar mass black hole. The key idea is 
that the relativistic outflow that we think is at the origin 
of the prompt gamma-ray emission, is nothing else 
than the magneto-centrifugally driven neutrino-wind, 
coming from the cooling proto-NS. The dynamics of 
such wind reaches rapidly the force-free limits, as the 
baryon loading rapidly drops in time (Metzger et al., 
20112), and its interaction with the surrounding layers of 
the progenitor star or the ejecta of the binary merger (for 
SGRBs), leads to the formation of a hot relativistic wind 
bubble, that is both a reservoir of energy, to be released at 
later time, as well as a reservoir of high energy photons 
that can lead to the appearance of a so called kilo-nova 
(Ciolfi & Siegel, 2015; Siegel & Ciolfi, 2016a,b). Lower 
energy systems have been considered also to justify 
a possible continuum of explosion phenomenology 
down to super-luminous and broad line Ib/c supernovae 
(Dessart et al., 2012; Moriya et al., 2017; Soker & Gilkis, 
2017; Milisavljevic et al., 2018; Wang et al., 2019; Vurm 
& Metzger, 2021; Shankar et al., 2021; Dessart, 2019; 
Margalit et al., 2018). 


3 MODELLING PWNE 


The pulsar wind of an oblique rotator has been shown to have 
an extremely complex structure named striped wind, with 
a magnetic field (B) of alternate polarities separated by a 
current sheet (Bogovalov, 1999). A very important parameter 
of the pulsar wind is the wind magnetization c, representing 
the ratio between the Poynting flux and the particle kinetic 
energy flux in the wind: 


B 
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with p the comoving density of the plasma and I, the Lorentz 
factor of the wind. In recent years our understanding of the 
properties of the pulsar wind has been increased thanks to 
very refined numerical models, going from force-free only 
(Spitkovsky, 2006, and subsequent) to more complex physics, 
including pair creation (Philippov & Spitkovsky, 2014). 

In this section we review the many different approaches 
that have been used to model PWNe, discussing merits and 
pitfalls of each approach, as well as their validity with respect 
to the different evolutionary phases. A complete description 
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of a PWN, to be compared with observations, requires the 
treatment of two different aspects: (1) the dynamical evolu- 
tion of the system; (2) the spectral evolution of the particles 
responsible for the emission. At present none of the models 
is able to account for both at a level of accuracy to enable a 
direct and reliable comparison with data. Multi-dimensional 
numerical simulations have reached in the last years excellent 
results in the description of the dynamics, but on the other 
hand they lack in terms of spectral modelling, with radiation 
generally evaluated on top of the dynamical results at the end 
of the computation, with simplified treatment of the radia- 
tion losses and of the particle energy evolution (Volpi et al., 
2008; Porth et al., 2014; Olmi et al., 2016). On the other hand, 
accurate radiative models only exist for very simplified 1D 
descriptions of the dynamics, thus the spectral evolution is 
pinned to a rough model of the physical properties of the 
system (Gelfand et al., 2009; Bucciantini et al., 2011; Torres 
et al., 2014). 


3.1 One-zone models 


In one-zone (or equivalently 0+1) models, the PWN is de- 
scribed as a uniform bubble in interaction with the surround- 
ing SNR, being subjected to adiabatic, radiation and possibly 
particles losses. Energy and magnetic field are injected into 
the bubble with constant relative efficiencies by the PSR, 
following its spin-down luminosity (Eq. 2). The magnetic 
field energy within the bubble can be evolved either assuming 
magnetic flux conservation (Pacini & Salvati, 1973; Gelfand 
et al., 2009; Torres et al., 2014) or a constant ratio between 
magnetic and plasma pressures (Bucciantini et al., 2011). In 
the early free-expansion phase, if radiation losses are negligi- 
ble, the two approaches are equivalent (however they might 
have different implications in terms of wind magnetization), 
but once radiation losses become important they might lead 
to sizeable different results (Bucciantini et al., 2011). 

One-zone models are all based on the thin-shell approxi- 
mation: the shell is evolved considering mass and momentum 
conservation (Eqs. 4-5 in Sec- 2.3) and its radius is equated 
with that of the PWN. To the time evolution of the radius, 
one then couples a description of the evolution of the inter- 
nal particles distribution function subject to injection and 
losses, from which one finally derives the observed PWN 
multi-wavelength spectrum. 

Spectral evolution models can be traced back to the original 
work by Pacini & Salvati (1973) that, assuming an approx- 
imated description of the temporal evolution of the PWN, 
limited to its initial free-expanding phase, derived the evolu- 
tion of the spectral energy distribution function showing its 
relation with the injected particle distribution. This work was 
later extended with a model also describing the interaction 
with the SNR, but still limited to the interaction with the 
unshocked ejecta (Reynolds & Chevalier, 1984), and consid- 
ering possible variations to the particles injection (Atoyan, 
1999). 

Despite their evident over-simplification, one-zone mod- 


els have been — and still are — widely used. Results from 
one-zone models have in fact proved to give a good descrip- 
tion of the global properties of young PWNe, allowing one 
to model the evolution of these systems and their spectral 
properties, and to sample a large parameter space, in ways 
that are not possible with current multi-dimensional models. 
They can be easy implemented and are not much demand- 
ing in terms of numerical resources (Gelfand et al., 2009; 
Tanaka & Takahara, 2010; Bucciantini et al., 2011; Martín 
et al., 2012; Tanaka & Takahara, 2013; Torres et al., 2014; 
Gelfand, 2017; Fiori et al., 2020). These characteristics make 
them very much appealing for large populations studies (e.g. 
Torres et al. 2014; Fiori et al. 2022), but on the other hand 
they must be used with care to follow the evolution beyond 
free-expansion. In Bandiera et al. (2020, 2023) it has been in 
fact shown that one-zone models predict excessive compres- 
sions of the nebula during reverberation, possibly leading to 
the burn-off of a huge amount of the emitting particles, chang- 
ing dramatically the spectral evolution. Torres et al. (2019) 
show that in cases of extreme compressions (CF z 1000) a 
super-efficient phase can appear, when the PWN luminosity 
is so enhanced to make it visible again at X-rays at later times. 
One-zone models tend to overestimate the PWN compression 
mainly because the simplified description of the pressure in 
the SNR, in general assumed to be equal to the central pres- 
sure from the Sedov solution, or the pressure at the SNR FS 
scaled with some arbitrary factor (for a complete description 
see Bandiera et al. 2023). This kind of approximation indeed 
introduces large errors when the PWN starts to interact with 
the SNR. Recently Bandiera et al. (2023) have in particular 
shown that the pressure in the SNR is very far from the Sedov 
solution during the entire first compression, and that in the 
aforementioned assumptions lead to a consistent overestima- 
tion of the outer pressure. This means that the number of 
systems undergoing a super-efficient phase, might be much 
smaller than expected, if a correct description of the SNR 
pressure is considered during reverberation. 

A preliminary attempt to extend one-zone models beyond 
free-expansion has been recently made by Fiori et al. (2022), 
where the SNR pressure is shaped on the results of 1D hy- 
drodynamic simulations of the PWN-SNR interaction. Con- 
versely to Bandiera et al. (2023), in this case there was no 
specific focus in reproducing correctly the dynamical effects 
of the interaction between the PWN and the surrounding 
SNR. 

The extension of one-zone models is generally limited 
to the first series of compressions and re-expansions of the 
PWN, while in few cases they have been extended to longer 
times halting by hands the oscillations to match the SNR 
pressure from the Sedov solution (Torres et al., 2019). 


3.2 1D models 


In 1D models the evolution and properties of the PWN and 
the SNR are given as a function of the radial coordinate. The 
first 1D model of a PWN was put forward by Rees & Gunn 
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(1974) to describe, as many others later on, the Crab nebula. 
Despite the very simplified description of the system, that 
model was able to predict a number of features and to give 
them a physical interpretation: the appearance of the ultra- 
relativistic un-shocked wind as an under-luminous area in the 
inner nebula, later observed by Weisskopf et al. (2000) at X- 
rays; the nebula shrinkage with increasing frequency as sign 
of the synchrotron emission and central particle injection; 
the nebular magnetic field close to the equipartition value; 
the Lorentz factor of the wind, its magnetization and the 
injection rate of particles from the synchrotron luminosity. 
This preliminary model was then elaborated and extended by 
Kennel & Coroniti (1984a,b), that provide estimates of many 
characteristic quantities of the Crab nebula based on the full 
solution of the relativistic MHD equations within the PWN 
itself. Few years later, Emmering & Chevalier (1987) found 
a time dependent analytic solution of the same problem. 

Many other 1D models have been then developed in later 
years, all based on the numerical implementation of the equa- 
tions describing the PWN-SNR interaction and evolution, 
both in the classic and relativistic regimes and in the hydro- 
dynamic (HD) or MHD frameworks, extending to a longer 
evolution than the free-expansion phase (van der Swaluw 
et al., 2001, 2004; Bucciantini et al., 2003; de Jager et al., 
2008; Bandiera et al., 2023). 

One of the longer standing problems in PWNe physics, 
the so called sigma-paradox (Melatos, 1998), arose as conse- 
quence of these 1D models, and puzzled the community for 
more than thirty years: in order to explain the existence of 
the TS, the average magnetization o in the nebula must be 
quite low: c = few x 10? (Kennel & Coroniti, 1984b). But 
theoretical models of pulsar magnetospheres predict a much 
larger magnetization at the light cylinder of the star: o ~ 10^ 
(Kirk et al., 2009; Arons, 2012). To make compatible these 
two opposite predictions, a huge amount of magnetic dissipa- 
tion must be considered along the path separating the light 
cylinder and the pulsar wind termination shock, converting 
the pulsar wind from a Poynting dominated outflow to a parti- 
cle dominated one. Magnetic dissipation is actually expected 
to occur in the current sheet of the alternating pulsar wind 
(Lyubarsky, 2003; Sironi & Spitkovsky, 2011), but ~ 7 orders 
of magnitude are simply too many to be accounted for with 
any known process. We want to remark here that the magneti- 
zation at the wind termination shock not only affects the PWN 
dynamics, but it is also an important parameter that controls 
the efficiency of the various acceleration mechanisms (Am- 
ato, 2014). As we will discuss in the following paragraphs, 
going multi-D is the way to reduce — if not solve - this long 
standing problem, since the augmented dimensionality allows 
increasing the amount of magnetic dissipation. 

Results of one-zone and 1D models are in excellent agree- 
ment in the free-expansion phase, while the accordance dis- 
appears as soon as the PWN enters reverberation. As for 
one-zone models, the reliability of 1D models is in any case 
limited to at most the beginning of reverberation. Unlike one- 
zone, however, 1D Lagrangian models (Bandiera et al., 2023) 
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can be extended to the following evolution (the first compres- 
sion and the subsequent sequence of multiple compressions 
and re-expansions). Recall that, in any case, these are just an 
artifact of the reduced dimensionality. 1D simulations show 
a sort of damped oscillations, with compressions generally 
becoming less severe as time passes by. These behaviour is 
not found in a multi-D description of the PWN/SNR system, 
where the mixing produced by boundary instabilities helps in 
mitigating oscillations. However, one should be careful not 
to confuse the observed angular size/dimension of the PWN 
with its effective volume (the one that matters for compres- 
sion), once mixing becomes important. 


3.3 2D models 


Overcoming the limitations of 1D models became urgent 
with the first detailed images of the inner Crab nebula at 
optical (Hubble Space Telescope, Hester et al. 1995) and 
X-rays (Chandra, Weisskopf et al. 2000). The Crab nebula 
was in fact the first source where a beautiful bright jet-torus 
structure, with an equatorial torus and two opposite polar jets, 
was identified. Later on other systems showed up the same 
properties (Gaensler et al., 2002; Lu et al., 2002; Romani & 
Ng, 2003; Camilo et al., 2004; Slane et al., 2004; Romani 
et al., 2005), now believed to be a common feature of young 
PWNe. It is clear that such a complex structure cannot be re- 
produced with the simplified geometry of 1D models, neither 
by slighting modifying the 1D approach. Indeed the first 2D 
analytic model of the Crab nebula by Begelman & Li (1992), 
considering the toroidal structure of the magnetic field, was 
only capable to explain its observed prolate shape. 

A critical point for the theoretical description of the inner 
nebula was the appearance of the polar jets so close to the 
pulsar. They in fact seem to form in the un-shocked relativis- 
tic wind, where magnetic collimation is known to be poorly 
efficient. Theoretical and numerical models of the relativistic 
wind emanating from the pulsar (Begelman & Li 1992; Con- 
topoulos et al. 1999; Bogovalov 1999, and later Komissarov 
2006; Spitkovsky 2006; Timokhin 2006) already predicted 
a non-uniform distribution of the wind energy flux at the 
termination shock, with most of the energy concentrated in 
the equatorial plane (the so called split-monopole models), 
but there was no evidence for a collimation of part of the flux 
in the polar regions. The solution to the jet formation was 
found only few years later: modelling the dynamics of the 
plasma with an anisotropic distribution of the energy flux in 
the wind, the terminations shock was found to become oblate, 
with larger extension in the equatorial region than in the polar 
one (Lyubarsky, 2002; Bogovalov & Khangoulian, 2002; Bo- 
govalov & Khangoulyan, 2002; Khangoulian & Bogovalov, 
2003). Polar jets can then form due to magnetic hoop stresses 
in the post-shock plasma, immediately beyond the polar front 
of the shock, and appear closer to the pulsar than the torus 
thanks to the oblate morphology of the TS itself. 

This theoretical prediction was verified later on with the 
use of relativistic 2D MHD numerical simulations (Del Zanna 
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et al., 2004; Komissarov & Lyubarsky, 2004; Bogovalov et al., 
2005; Del Zanna et al., 2006). Moreover, 2D models show 
that the formation of jets is only possible starting from a 
minimum wind magnetization in the pulsar wind of o ~ 1072, 
one order of magnitude larger than that originally set by 1D 
models. The increase of dimensionality then appears as the 
first possible way to alleviate the sigma-paradox. 

Thanks to its luminosity and vicinity, the largest part of 2D 
models were made to investigate the properties of the Crab 
nebula that was — and still is — the perfect source to look at 
for a detailed comparison. They were extremely successful 
at reproducing many of its features down to very fine details, 
especially in the inner nebula: maps of the synchrotron emis- 
sion (Del Zanna et al., 2006; Olmi et al., 2014); polarization 
properties (Bucciantini et al., 2005); the complete spectrum, 
from radio to gamma-rays (Volpi et al., 2008); the variability 
at small scales in the inner nebula (Camus et al., 2009; Olmi 
et al., 2015; Lyutikov et al., 2016). In particular the multi- 
wavelength appearance of variable arc-like bright structures 
(the so called wisps), forming close to the TS and moving out- 
wards at a consistent fraction of c, was shown to be a perfect 
tool to trace the properties of the underlying plasma and to 
constrain the location and mechanism of particle acceleration 
(Olmi et al., 2015). 

On the contrary, results of 2D models at larger scales are 
limited by the imposed axisymmetry, with only the global 
properties of the nebula, such as the extension and the shrink- 
age with the increasing energy, reproduced reasonably well. 
Axisymmetry is in fact the major limitation of 2D models: it 
induces the artificial accumulation of magnetic loops along 
the polar axis of the nebula, producing enhanced compression 
of the magnetic field. This imposes an upper limit to the wind 
magnetization of c < 0.1, otherwise the polar compression of 
the field becomes so strong that the collimated magnetic flux 
punches the nebular shell at the polar boundaries. The con- 
sequence of this non-redistribution of magnetic field is that 
polar regions have an (excessively) intense magnetic field, 
while the rest of the nebula is under-magnetized, with an 
average magnetic field, for the Crab, well below the expected 
value of ~ 150 uG. To reproduce the PWN spectrum then one 
is forced on one side to inject an excessive number of parti- 
cles, on the other to steepen by hand the injection spectrum of 
the X-ray emitting component to alleviate the average lower 
energy losses, and match the synchrotron emission. 

The problem introduced by the wrong geometry of the 
magnetic field remains very evident if looking at the mor- 
phology of the nebula at large scales: the radio emission, 
expected to be rather uniform, indeed resembles the shape of 
the magnetic field (it is concentrated around the polar axis); 
the IC emission is overestimated, approximately by the same 
factor by which the averaged magnetic field is underestimated 
(Olmi et al., 2014). 

Despite these limitations, 2D models remain appealing in 
many cases: their results are in good agreement with more 
complex 3D models if limited to the inner nebula and small 
scales; they allow for a longer evolution of the PWN and for 
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a larger investigation in terms of number of sources and phys- 
ical parameters with respect to their equivalent in 3D, due 
to the lower numerical cost. A possibility to use 2D models 
to infer the large scale properties of PWNe has been inves- 
tigated in Olmi & Torres (2020), where the HD scheme has 
been preferred to the MHD one to avoid the aforementioned 
problems with the field geometry, affecting the large scale 
emission. In that work the magnetic field is excluded from 
the dynamics but traced numerically (with a recipe linking 
it to the thermal pressure), approximately accounting for the 
particles losses during the evolution. 

Finally, 2D HD and MHD models have also been suc- 
cessfully used to investigate the formation of bow shock 
PWNe produced by fast moving pulsars escaped from their 
SNRs, and their variations depending on the properties of the 
ambient medium (Toropina et al., 2001; Bucciantini, 2002; 
Bucciantini et al., 2005; Olmi et al., 2018; Toropina et al., 
2019). 


3.4 3D models 


With the first 3D MHD simulations of the Crab nebula (Porth 
et al., 2013) it became apparent that the solution to the sigma- 
paradox was finally in reach. We have already seen how 
moving from 1D to 2D permitted to gain more than one 
order of magnitude in wind magnetization. Thanks to the 
development of efficient processes of magnetic dissipation in 
the nebula, a magnetization c > 1 becomes finally accessible 
with 3D simulations, drastically reducing the impact of the 
sigma-paradox. Configurations with a toroidal magnetic field 
are subject to current driven instabilities (kink-like, Begelman 
1998; Nalewajko & Begelman 2012), and signs of such a 
process have been detected in different PWNe (Mori et al., 
2004; Pavlov et al., 2003). Numerical simulations of the Crab 
nebula jet proved not only that the kink instability efficiently 
develops in the nebula, but that it is also responsible for 
the variations seen in the jets at different epochs (Mignone 
et al., 2013). Mixing of the magnetic field induced by kink- 
like instabilities is so efficient that, even if the initial field 
configuration is fully toroidal, a poloidal component raises 
rapidly in 3D, becoming even dominant immediately outside 
the inner nebula (see e.g. Porth et al. 2014; Olmi et al. 2016). 
The mixing causes on one side the magnetic field geometry 
to become much more complex than in 2D, but also results 
in a more uniform distribution of the magnetic field in the 
nebular volume (see Fig. 4). 

Despite being the best way to account for the properties of 
a PWN, the use of 3D models is limited by the huge amount 
of resources (time/numerical) they need. The spatial scales 
that must be reproduced are extremely different, from the 
injection region of the pulsar wind (smaller than the TS) to 
the PWN contact discontinuity, and the surrounding SNR 
bubble in case of young systems. Moreover to correctly repro- 
duce the pulsar wind, the injection site must be solved with a 
minimum number of grid cells (z 10 — 20), and the injection 
region must always remain detached from the radius of the 
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Figure 4. 3D plot of the magnetic field lines (left side) and slice of the magnetic field intensity (right side) from the 3D simulation of the Crab nebula presented 
in Olmi et al. (2016), at a simulated age of ~ 250 yr. Stream lines are colored with the intensity of the field as from the right panel, but using a different color 
scale to highlight the different structures. Both figures have been produced using VisIt (Childs et al., 2012). 


termination shock, otherwise the correct jump conditions at 
the shock cannot be ensured. This translates in the necessity 
of a very high resolution at the center of the numerical do- 
main, mapping a zone which represents only ~ 1/100 of the 
global size of the system. The grid is then usually optimized 
with the use of an adaptive mesh technique, able to increase 
or decrease the resolution as needed, allowing to save a large 
amount of time, or with expanding grids. Nevertheless the 
resources needed to run such models remain large and still 
prohibitive for long term evolution. To date the longest 3D 
MHD simulation of the Crab nebula reproduces ~ 1/4 of the 
age of the source (Olmi et al., 2016) and it required few mil- 
lions of core/hours of computational time and several months 
to be run. A longer evolution (~ 7500 yr) has been reached 
with HD simulations and expanding grids in the case of a high 
speed PWN (vxick ~ 300 km s^!) interacting with a composite 
SNR (Kolb et al., 2017), producing a strongly asymmetric 
system out of the reverberation phase. 


A different approach has been applied to the modelling of 
bow shock pulsar wind nebulae. For these systems, going 3D 
is particularly important to correctly capture the magnetic 
field topology and the development of turbulence in the tail, 
which are strongly affected by geometric constraints and 
relevant for the observational properties of the source. A 
first attempt to 3D modelling of bow shocks, limited to the 
classical HD regime, was presented by Vigelius et al. (2007). 
Then, in recent years, the extension to the MHD relativistic 
regime was investigated by different groups (Barkov et al., 
2019; Olmi & Bucciantini, 20192). In all these cases, the bow 
shocks have been modelled directly assuming the pulsar in 
interaction with the ISM, not considering then the transitional 
phase when the star is emerging from the SNR bubble. That 
transitional phase was indeed investigated in van der Swaluw 


et al. (2003). 

The PSR reference frame is generally used, with the star 
being positioned at a specific point of the computational 
domain. The star then sees the ambient medium as an incom- 
ing, cold flow, that might be magnetized or un-magnetized, 
depending in the specific case. A large sample of different 
configurations has been investigated, varying the inclination 
of the magnetic field and pulsar-spin axis, the direction of the 
pulsar motion compared to the first two, the magnetization 
and level of anisotropy of the pulsar wind. The evolution is 
then followed for enough time to ensure the system dynamics 
has reached a relaxed quasi-steady regime. Emitting proper- 
ties can be computed using the same approach generally used 
for young systems. Particles responsible for the emission are 
injected at the wind termination shock with a broken power 
law distribution and their emissivity is computed as discussed 
in the following Sec. 4, with different possible assumptions 
on their density in the PWN (e.g. they can be considered as 
uniformly distributed or with a distribution shaped on the 
thermal pressure, see Olmi & Bucciantini 2019b). 

In a recent work (Olmi & Bucciantini, 2019c) we have 
addressed the possibility for particles to escape the bow shock, 
modelling the evolution of the particle trajectories in the 
electric and magnetic fields of the MHD simulations. 


3.5 Highlights from 3D MHD models of young PWNe 


Here we summarize the most important findings of 3D nu- 
merical simulations of young — Crab like - PWNe, according 
to Porth et al. (2013, 2014) and Olmi et al. (2016). 


e Solution of the sigma-paradox: values of the wind mag- 
netization at injection larger than unity become possible 
thanks to the efficient magnetic dissipation produced in 
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Figure 5. A selection of images for a PWN in the free-expansion phase, with maps from MHD simulations specific for the Crab nebula, elaborated from the 
simulations described in Del Zanna et al. 2006 (middle panel — 2D) and Olmi et al. 2016 (panel on the right — 3D). In particular, from left to right: toy model 
of the inner structure of the nebula, showing the oblate shape of the wind termination shock (the yellow spiral marks the striped wind) and the formation of jets 
due to hoop stresses at the polar front of the shock; 1 keV X-ray synthetic map of the Crab nebula, reproducing its inner morphology (map in logarithmic scale 
and intensity scaled to the maximum); intensity map of the 3D velocity field from a 3D simulation of the Crab nebula, with the evident formation of kinking 
jets. 


3D. Actually the average value of the magnetic field when 
the system has reached a self-similar evolution (at ages 
= 150 yrs) is still a factor of ~ 1.5 — 2 lower than the 
expected one (see Fig. 6 in Olmi et al. 2016). 
Development of poloidal magnetic field: despite the mag- 
netic field at injection is purely toroidal, turbulence and 
high-speed polar flows rapidly modify its topology. A 
polar component easily develops immediately beyond the 
inner nebula, and becomes even dominant in the outer 
nebula (see e.g. Fig. 11 in Porth et al. 2014 or Fig. 8 in 
Olmi et al. 2016). The complex structure of the magnetic 
field topology and its intensity can be seen in Fig. 4. 


found in the outer nebula, very similar to what observed 
at large scales in radio (Bietenholz et al., 2004). 

Lacking of correct spectral description: despite the huge 
improvements in the description of the fluid structure, 
and in general in the global dynamics of the nebula, the 
emitting properties are not fully reproduced, both in terms 
of morphology and integrated spectrum. The cause may 
be both in the — still — too low magnetic field on average, 
and/or in the approximated reconstruction — a posteriori — 
of the properties and history of the emitting particles. 


e Despite the complex topology of the magnetic field, the 3.6 Highlights from 3D MHD models of bow shocks 


magnetic pressure in the nebula is rather uniform on large 
scales (see Fig. 9 in Olmi et al. 2016), as well as the to- 
tal pressure (magnetic and thermal, see Fig. 3 in Porth 
et al. 2014). This is a net difference with 2D MHD ax- 
isymmetric models, and explains why an HD approach 
seems more suitable to reproduce bulk and macroscopic 


In this subsection we summarize the results of 3D relativistic 
MHD simulations of bow shock nebulae. The discussion is 
mostly based on Barkov et al. (2019) and Olmi & Bucciantini 
(2019a) for the dynamics, on Olmi & Bucciantini (2019b) 
for the emission and polarization properties and on Olmi & 
Bucciantini (2019c) for the properties of the particle escape. 


properties of PWNe (Olmi & Torres, 2020). 

The inner nebula is reasonably well described by 2D 
MHD axisymmetric models, including the Crab wisps 
and knot (Komissarov & Lyubarsky, 2004; Lyutikov et al., 
2016). A direct comparison shows that, if one sticks to the 
inner — toroidal — nebula, the description obtained with 2D 
models does not differ much from the structure simulated 
in 3D (see Fig. 5). This of course is a very important 
result, supporting the reliability of previous models in 2D 
and meaning that, if only interested in the properties of 
the inner nebula, one can safely use less demanding 2D 
simulations. 

Time variability of the inner nebula: the wisp-like vari- 
ability close to the shock position is reproduced correctly 
and is comparable with what previously found with 2D 
models (Camus et al., 2009; Olmi et al., 2015). In 3D a 
number of variable non-axisymmetric structures is also 


e The large scale structure of the bow shock is rather inde- 


pendent on the variation of the geometry and wind proper- 
ties, with major deviations caused by different inclinations 
of the pulsar spin-axis and the direction of motion (maxi- 
mum for 45°). The differences appear as small extrusions 
and blobs, in some cases resulting in periodic perturbation 
of the forward shock. 

The dynamics in case of an isotropic wind is very similar 
to the one found in 2D HD for intermediate values of the 
magnetization (o = 0.1). 

The dynamics in the tail is largely dominated by the 
level of magnetization and by the wind anisotropy: unlike 
anisotropy, the lower is the magnetization, the higher is 
the turbulence. Anisotropic models then are more turbu- 
lent than isotropic ones, showing strong mixing also for 
large magnetization (see Fig. 3 of Olmi & Bucciantini 
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Figure 6. Comparison of two different regimes of particle escaping from a 
bow shock nebula from a 3D simulation: from massive and almost diffusive 
escape for high Lorentz factor particles (y = 108, left-blue colored side) 
to directional escape along external field lines of lower energy particles 
(with y = 107, right-orange colored side). Data come from the simulation 
presented in Olmi & Bucciantini (2019c) and have been displayed using 
VisIt (Childs et al., 2012). The inclination of the image is shown with the 
bottom triad, while the position of the bow shock is marked by the dashed 
white line. 


20192). In general, anisotropic and low magnetized sys- 
tems are fully dominated by turbulence even close to the 
termination shock, resulting in the complete loss of in- 
formation from the injection region along the tail. On 
the other hand, isotropic and highly magnetized systems 
show a coherent structure of the magnetic field, with the 
injection properties still affecting the dynamics far from 
the shock along the tail, in very good agreement with 
simplified semi-analytic models (Bucciantini, 2018). 

e Noefficient magnetic amplification from turbulence. Even 
in presence of strong turbulence, the magnetic field tends 
to reach equipartition with the turbulent kinetic energy, 
in general smaller then the thermal energy in the tail. 
The only sign of field enhancement is found close to the 
contact discontinuity, possibly resulting from efficient 
shear instability amplification. 

e Low turbulence cases — emission and polarization proper- 
ties: a strong correlation between conditions at injection 
and surface brightness is found. The variety of observa- 
tional morphologies is very wide: from bright heads to 
bright tails or, in some cases, bright wings. The polariza- 
tion fraction is higher in the tail for higher magnetization. 

e High turbulence cases — emission and polarization prop- 
erties: once magnetization drops (and anisotropy grows), 
turbulence starts to dominate also in the appearance of 
surface brightness, and distinguishing between the vari- 
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ous cases becomes hard. When turbulence increases, the 
polarization fraction drops. 

e The escape of particles from the bow shock is found to be 
an energy dependent process: the threshold for escape is 
set by the condition that the particle Larmor radius (rz) in 
the equipartition magnetic field (Beg ~ few uG) is equal 
to the typical size of the bow-shock in the head (the stand 
off distance do). This translates in a Lorentz factor of the 
particles of: y ~ eBeg do/ (m,c?), that for typical systems? 
corresponds to y = few x 10’, or in a particle energy = 10 
TeV. 

e There is a transition in the escape process: particles man- 
age to escape more easily if injected at the frontal polar 
region of the pulsar wind, while the others tend to remain 
confined in the tail. At lower energies, particles escape 
only in the presence of reconnection points at the magne- 
topause between the shocked pulsar wind and the ISM, 
and this might give rise to the appearance of one-sided 
jet-like features. Particles show an increasingly more dif- 
fusive escape with energy: the outflow becomes more 
uniform (see Fig. 6), but charge separation increases. 


4 RADIATION AND ACCELERATION 


The pulsar is an excellent conductor. Charges inside it orga- 
nize themselves in such a way that the internal electric field 
is screened. But the electric field at the star surface is not, and 
it is strong enough to extract charged particles (leptons and 
possibly ions) from the star. This generates a co-rotating mag- 
netosphere that extends up to the star light cylinder Rzc = cP, 
with P the pulsar period. The magnetic field lines originating 
close to the pulsar magnetic axis (at the so called polar caps) 
extend beyond the light cylinder and form the open magneto- 
sphere, through which the pulsar wind flows into the nebula. 
The rate at which particles are extracted at the PSR surface is 
given by Goldreich & Julian (1969): 
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with By. the magnetic field at the polar cap and ® = VE/c 
the maximum potential drop between the pulsar and infinity. 

Once extracted, particles are accelerated at different loca- 
tions along the open field lines, where they meet regions of 
un-screened electric potential. As a consequence they emit 
high-energy photons that can be absorbed in the intense mag- 
netic field surroundings the star, to generate electromagnetic 
cascades. The number of pairs in the magnetosphere then 
increases by a large factor measured by the pair multiplicity 
k ~ 10^ — 107, namely the number of secondary leptons gen- 
erated from the primary extracted from the star. The exact 
estimate of the multiplicity is very controversial (Timokhin 
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in units of 10°° erg s~!, the ambient density in units of 1 proton per cm? and 
the PSR velocity in units of 200 km s^. 


The evolutionary path of PWNe 


& Harding, 2019), and there is a possibility to infer it from 
the modelling of the PWN properties. 

Differently from leptons, ions cannot be generated in cas- 
cades, and so, if present, they must be a factor of k less than 
pairs. But given the difference in mass between ions and elec- 
trons (m,/m, ~ 1800), this does not necessary means that 
they are irrelevant in the PWN energetics (see e.g. Amato & 
Olmi 2021 and the discussion therein). 

PWNe reprocess a consistent part of the pulsar spin-down 
power into accelerated particles, with the Crab being the most 
efficient known at ~ 30% efficiency. Only a very small frac- 
tion goes into pulsed radio to X-ray radiation (x 1%), while 
a larger one might go into pulsed gamma-rays (Abdo et al., 
2013a). In general the study of the PWN emission is then 
relevant to obtain indirect information about pulsar physics. 
PWNe shine at multi-wavelengths via non thermal emission. 
The primary emission mechanism, at least for a consistent 
period of their life (in free-expansion and possibly for large 
part of reverberation), is synchrotron radiation produced by 
the shocked wind particles interacting with the nebular, rather 
intense (^ 50 — 200 uG), magnetic field. The synchrotron 
spectrum can be modelled as a set of broken power-laws. 
One (or two, for t >> To, see Pacini & Salvati 1973) break is 
associated to synchrotron cooling (E<). However, the others 
are typically thought to be associated with changes in the 
acceleration mechanism. The exact location of these breaks 
in the spectrum cannot be trivially inferred, and it is often 
based on more detailed modelling. In the Crab nebula, current 
models suggest that the one between optical and X-rays is 
due to synchrotron cooling, while the one between radio and 
optical is attributed to a change in the particles acceleration 
mechanism. In other systems like MSH 15-52 the two breaks 
are so close, and the spectral coverage so sparse, that it is 
hard to guess what is what (Nakamori et al., 2008; Gaensler 
et al., 1999, 2002). 


4.1 Injection at the shock 


The particle injection spectrum is thus typically modeled as a 
broken power-law in the particle energy E: 
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consisting of two distinct families, characterized by different 
injection indices p;, at energies below or above the injection 
break (Ep). The normalization function Qo(t) is determined by 
the requirement that the power injected in particles is a fixed 
fraction of the spin-down luminosity, namely: (1 — 7)L() = 

p Q(E,t)E dE, where 7 is the magnetic fraction (i.e. how 
much of the injection goes into magnetic energy), linked to 
the magnetization (Eq. 13) as: 7 = c/(o + 1). 

A rather flat injection index (pios ~ 1 — 1.5) is characteris- 
tic of the lower energy component, leading to a synchrotron 
spectrum with flux density S ,(v) x v *' and spectral index 
Qoy ~ 0 — 0.3. The higher energy part indeed shows steeper 
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spectra at injection (pio, ~ 2 — 2.7), leading to a synchrotron 
spectrum with @pjgn ~ 1 — 1.2 or, as more commonly reported, 
a photon index I = 1 + œ ~ 2 — 2.2. The recent PeV observa- 
tions by LHAASO (Cao et al., 2021a,b) confirm, as originally 
suggested by Bucciantini et al. (2011) and recently inves- 
tigated by de Ofía Wilhelmi et al. (2022), that high energy 
particles can be accelerated up to the pulsar voltage (see also 
Khangulyan et al. 2020). 

The possibility that particles responsible for the radio, op- 
tical and X-ray emissions belong to separate families, ac- 
celerated via different mechanisms, is supported also by the 
multi-wavelength variability observed in the Crab nebula. The 
arc-like bright structures named wisps, that appear very close 
to the TS location and move outwards with mild relativistic 
velocity, have been observed at multi-wavelengths (Hester 
et al., 2002; Bietenholz et al., 2001, 2004), and shown to 
neither be spatially coincident nor characterized by the same 
velocity (Schweizer et al., 2013). In the MHD framework, 
which has provided a very good description of the nebular 
dynamics, wisps trace the structure of the underlying plasma 
— the magnetic field in particular — and as such they can only 
be non-coincident if particles with different energies are pro- 
duced at different locations of the shock. This likely means 
that acceleration processes act differently in different sectors 
of the shock (Olmi et al., 2015). 

The acceleration mechanisms proposed so far are mainly 
three: (1) diffusive shock acceleration, or Fermi-I like pro- 
cesses; (ii) diffuse acceleration due to stochastic magnetic 
reconnection, or Fermi-II, in MHD turbulence; (iii) accelera- 
tion conveyed by driven magnetic reconnection. Actually a 
fourth mechanism has been invoked: resonant absorption of 
ion cyclotron waves, which however requires the presence of 
ions in the wind, still not confirmed (nor excluded). 

Diffusive shock acceleration requires a very low magne- 
tization to be effective (c € 10-3, Sironi et al. 2015), a 
condition that can only be sustained at the equatorial sector 
of the shock, where the striped wind ensures a huge dissipa- 
tion of the field, or very close to the polar axis, where the 
field naturally vanishes. The power law index at injection 
of optical/X-ray emitting particles is compatible with what 
predicted by Fermi-I acceleration, and MHD models for the 
X-ray wisps are also in agreement with a scenario in which 
that particles are injected mainly at the equatorial front of 
the oblique termination shock. On the other hand, driven 
magnetic reconnection requires a much higher magnetization 
(c z 30) and a large pair multiplicity («x = 10? ), difficult 
to account within present models of pulsars magnetospheres 
(Sironi & Spitkovsky, 2011; Timokhin & Harding, 2019). The 
power law index of the radio emitting particles at injection 
is on the other hand compatible with both reconnection and 
Fermi-II acceleration, while no particular information arises 
from wisps models in this case. A possible conclusion is that 
radio particles are accelerated at higher latitudes along the 
shock, out of the equatorial sector, or via a more distributed 
acceleration in the nebula (Olmi et al., 2015; Lyutikov et al., 
2019). 
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4.2 Radiation mechanisms 


Once injected in the nebula, the particle energy E evolves 
according to the following equation: 
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where c2 = 4/3 Oom/ (mc), with c the Thomson cross sec- 
tion, and U;aq is the energy density in radiation. The first term 
on the right side of the equation represents the energy vari- 
ation due to adiabatic expansion or contraction (depending 
on the evolutionary phase), while the second one describes 
radiation losses (synchrotron and IC). From this equation one 
can easily find that the cooling energy Ee is given in general 
by: 

1 ôR 1 

R Ot c3 (B2/8z + Uraa) ` 
Particles with E < E. are then dominated by adiabatic pro- 
cesses, losses or gains depending on the expansion or com- 
pression of the PWN. All particles with energy above E, are 
instead dominated by radiation losses. In the free-expansion 
phase E, determines the cooling spectral break, that it is com- 
monly found to move to higher energies with time (e.g. Pacini 
& Salvati, 1973). During reverberation instead E, separates 
those particles that gain energy due to compression (E < E.) 
from those that are still loosing energy due to radiation losses 
(E > Eo). 

As already mentioned, the spectral energy distribution of a 
PWN is fully non-thermal, with the primary emitting mecha- 
nism, responsible for emission from radio up to few hundreds 
of MeV, being synchrotron radiation. Higher energies are 
produced with the second emitting process characteristic of 
those systems: IC scattering between local photons and the 
same leptons responsible for the synchrotron emission. An 
example of a full spectral energy distribution coming from a 
one-zone modelling can be found, in the specific case of the 
Crab nebula, in Fig. 7. 

The main contribution to IC in general comes from the 
interaction with the photons of the cosmic microwave back- 
ground (CMB), with minor contributions from the interstellar 
radiation field and the synchrotron photons from the PWN 
itself. An exception in this respect is the Crab nebula: due to 
its young age and intense magnetic field (~ 150 — 200 uG), 
its IC spectrum has in fact a consistent contribution from self- 
synchrotron radiation. The Crab nebula is moreover the only 
known case where the maximum energy of the accelerated 
particles is limited by radiation losses in the Galaxy, and is 
smaller than the maximum energy inferred by the available 
potential from the pulsar. 

Once the injection spectrum has been defined (eq. 15), the 
PWN luminosity can be computed as: 


f ooa, a 


where Vpww is the PWN volume, while j? is the emis- 
sivity, obtained integrating over the particle distribution 
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Figure 7. Non-thermal spectral energy distribution of the Crab nebula, elab- 
orated from the original plot in Bucciantini et al. (2011), computed with a 
one-zone radiative model. Details of the assumed parameters and origin of 
data points can be found in the reference paper (see Fig. 1 and its caption). 
The different colored areas highlight the emission at various energy bands. 
From left to right: light-red for radio and IR; light green for optical and UV; 
light blue for X-rays (considering 0.1-100 keV); light yellow for the low 
energy gamma rays (up to the synchrotron limit ~ 250 MeV) and darker 
yellow for the high energy gamma-rays (fully due to IC). The blue line is 
for the synchrotron component, the red one for the IC component, to which 
major contributions come from self-synchrotron Compton (in cyan) and 
scattering with CMB photons (in yellow). 


function either for synchrotron or IC radiation as p © = 
T ™ OCE. t) P(E, y) dE, where Õ(E, t) is the evolved par- 
ics spectrum in the nebula, from the injection one defined 
in Eq. 15 (see e.g. Bucciantini et al. 2011). When computing 
the emissivity, in the 3D case, one should of course take also 
into account the spatial dependence due to orientation of the 
line of sight. General expressions for the synchrotron and IC 
power P(E, v) can be found in many textbooks, e.g. Rybicki 
& Lightman (1979). 


5 OBSERVING PWNE 


The firmly identified PWNe, with a detected associated pul- 
sar, to date counts ~ 60 systems. Most of them have been 
detected at X-rays thanks to Chandra (Kargaltsev & Pavlov, 
2008), while ~ 30 are those detected also at gamma-rays with 
different instruments (see e.g. Kargaltsev et al. 2013 and the 
TeV Cat catalog? for an updated list). The number of PWNe 
detected up to now in the Galaxy can be as high as ~ 90 
if we also include those sources marked as putative PWNe 
from their spectral and morphological properties, but with no 
associated pulsar. Another factor of ~ 4 is gained if we con- 
sider that a large part — if not all — of the present unidentified 
sources in gamma-ray surveys (e.g. Abdalla et al. 20182) are 
believed to be PWNe that have not been detected at lower en- 
ergies, possibly due to their evolved stage and the consequent 
faint emission at lower energies. If we consider a rate of birth 


The latest version of the TeVCat catalog con be found here: tevcat2. 
uchicago.edu 
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Figure 8. Distribution of the pulsars associated with identified PWNe on 
top of the complete pulsar population (in light blue) as taken from the 
ATNF catalog (Manchester et al., 2005), version 1.67. X-ray detected PWNe 
are shown as cyan circles or blue circles, in the last case those associated 
with fast moving pulsars. The Crab pulsar is shown as a red circle. All 
systems are given, respectively, in Table 1 and 2 of Appendix A, with some 
useful parameters. For an easier interpretation of the plot we also give lines 
indicating the range of the surface magnetic field characteristic of pulsars 
associated with PWNe (in light blue, Kargaltsev & Pavlov 2008), in gray 
lines of characteristic age from 10? yr to 10? yr and in pink lines of fixed 
spin-down luminosity 105? — 1036 — 109? erg s-!. 


of 107? pulsars yr^! in the Galaxy (Faucher-Giguére & Kaspi, 
2006), and an estimated lifetime of 10? yr at gamma-rays, the 
expected number of detectable PWNe at these energies in the 
Milky Way is in fact much larger than those really observed, 
with around ~ 1000 the expected systems. 

In Fig. 8 we mark the position in the P — P diagram of 
the known PWNe (i.e. of their associated pulsar), to be com- 
pared with the total population of pulsars as derived from 
the ATNF catalogue (Manchester et al., 2005, version 1.67, 
counting around 3300 stars). It can be noticed that PWNe 
appear to be associated only with the youngest part of the 
pulsar population, with age ranging between a few hundreds 
of years to a million of years. This can of course be partially 
due to a bias introduced by our inability to properly identify 
evolved systems, as in the case of the PWNe hidden in the 
population of unidentified gamma-ray sources, lacking of a 
multi-wavelength association. The lifetime of a PWNe as an 
X-ray synchrotron nebula is in fact much smaller than its life- 
time at very high energies. A 1 TeV photon can be produced 
via IC from the CMB - or possibly from the IR background 
— by an incoming electron of ~ 10 TeV energy, while the 
same photon requires a much more energetic particle to be 
produced via synchrotron radiation in the typical magnetic 
field of a PWN. In fact, even considering a very low magnetic 
field of ~ 10uG, characteristic of evolved systems, a 50 TeV 
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electron is necessary to produce a 1 keV photon. 

The lifetime of a lepton of energy Ee,Tey, expressed in units 
of TeV, against synchrotron losses in a magnetic field Bc, in 
units of uG, is in fact given by: 


-2 zj 
~ 25 EJ E yr. (19) 


It is then clear that the more energetic the lepton is, the 
quicker it radiates its energy away in the form of synchrotron 
emission, and the less long it survives. The same can be seen 
if looking instead at the energy of the synchrotron emitted 
photons (in keV units): 


Bao 7 { Egnxev V2 
Tsynch = 55.2 (5) (==) yr. Q0) 
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This makes a PWN detectable at X-rays only for a limited 
fraction of its life, when the pulsar is still powerful enough. 
On the contrary a PWN shines at radio energies for longer 
time, being the lifetime of radio emitting electrons much 
longer. These are also the same particles responsible for the 
long-living IC gamma-ray emission. As discussed previously, 
the detection at radio frequencies, especially for evolved, 
extended or diffused systems, might be difficult for multiple 
reasons, first of all instrumental limitations. Old nebulae are 
then likely to be detected mainly at gamma-rays, where we 
still lack in resolution, and their morphology is then difficult 
to be determined. 

A high level of linear polarization is one of the key proper- 
ties of synchrotron emission (Westfold, 1959; Legg & West- 
fold, 1968). For the typical particles distribution functions 
that are observed in PWNe, the polarized fraction theoreti- 
cally can be as high as 70%. It was indeed thanks to its high 
optical polarization, that synchrotron emission was recog- 
nized for the first time as the main emission mechanism in an 
astrophysical source, the Crab nebula (Baade, 1956; Oort & 
Walraven, 1956; Woltjer, 1958; Velusamy, 1985). 

Polarization is customarily observed in radio, and maps are 
available for many PWNe. The naive expectation is that the 
polarized structure in PWNe should correspond to a mostly 
toroidal magnetic field, as the one generated by a fast spinning 
rotator. There are indeed a few systems like Vela (Dodson 
et al., 2003) and G106.6+29 (Kothes et al., 2006a) where 
a well defined large scale toroidal pattern is observed with 
polarized fraction as high as 30%-40%. However, there is a 
wide variety in the radio polarization structures: some sys- 
tems show a large scale radial/dipolar pattern (Kothes et al., 
2008; Lai et al., 2022); while others have more random one, 
like the Crab (Bietenholz & Kronberg, 1990; Aumont et al., 
2010), with little to no correlation with respect to bright emis- 
sion features. Polarization is also available for old systems 
(Ma et al., 2016) and for a handful of bow-shocks (Ng et al., 
2010; Ma et al., 2016). Being radio emission in young or mid- 
dle aged PWNe dominated by the outer regions (since radio 
emitting particles are older and then fill the entire nebula), 
more subject to the interaction with the environment, radio 
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polarized measures provide at best a good estimate of the 
degree of ordered versus disordered magnetic field for the 
overall nebula, but cannot be used to investigate the condi- 
tions in the inner regions, where particle acceleration takes 
place. 

Optical and near IR polarization is only available for three 
systems: the Crab (Hester, 2008; Moran et al., 2013) where, 
due to the presence of a large foreground, only the bright 
knot and wisps have been studied and show a high level of 
polarization of 4096-5096, compatible with a toroidal mag- 
netic field; G21.5-0.9 (Zajczyk et al., 2012), where a small 
internal torus is observed with polarized fraction as high as 
50%; SNR 0540-69 (Lundqvist et al., 2011). 

Until the launch of the Imaging X-ray Polarimetry Explorer 
in December 2021 (IXPE, Weisskopf et al., 2022), the Crab 
nebula was the only PWNe (in-fact the only astrophysical ob- 
ject) to have a measured X-ray polarization (Weisskopf et al., 
1978). The polarized fraction was found to be 19%, with a po- 
larized angle marginally compatible with the symmetry axis 
inferred from fitting the X-ray torus (Ng & Romani, 2004). 
With IXPE an X-ray spatially resolved polarized measure 
finally became available for Crab (Bucciantini et al., 2022), 
Vela, and MSH 15-52, while integrated polarimetry will also 
be measured in a handful of other PWNe. 

In Appendix A an updated "catalog" of all the Galactic 
PWNe with an associated PSR known at present, divided in 
low and high speed systems, is reported in the two tables, 
with some ancillary information. 


5.1 Young systems 


To date we have identified less than 20 sources as PWNe in 
their free-expansion phase. They constitute a large part of 
the catalog of the PWNe detected mainly at X-rays and not 
associated with a fast moving pulsar (~ 40 sources), reported 
in Table 1. From the morphological point of view, synchrotron 
dominated systems are characterized by a larger extension 
at lower energies than at higher ones (see e.g. Fig. 9), with 
sort of a spherical/elliptical shape. X-rays highlight the inner 
nebula, revealing the presence of a jet-torus structure (see the 
zoom in of Fig. 9), believed to be a rather common feature in 
young PWNe, first discovered in the Crab and successively 
identified in another bunch of systems. 

Another common feature of the observed torii is the ap- 
pearance of enhanced brightness at one side (effect of the 
relativistic Doppler boosting of the emitting particles moving 
towards the observer) and the presence of variable, both in 
brightness and position, arc-like (the wisps, Scargle 1969; 
Bietenholz et al. 1991, 2001; Helfand et al. 2001; Pavlov et al. 
2001; Bietenholz et al. 2004) or point like (knots, Lou 1998) 
structures marking the high variability of the inner nebula, 
where (most of?) the particles are accelerated. As mentioned 
in section 3.3, the discovery of the complex inner structure 
of the Crab nebula was what prompted the move from 1D 
models to 2D MHD simulations. 

As discussed previously (Sec. 4), young PWNe are char- 
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acterized by extremely broad band spectra, extending from 
radio to gamma-rays; with the advent of LHAASO, now the 
Crab spectrum has been further extended above PeV energies 
(Cao et al., 202 1a). 


5.2 Middle aged - reverberating systems 


The reverberation phase, characterized by the interaction with 
the SNR reverse shock, is hardly identifiable. At present we 
only know a handful of systems showing clear evidence of 
being in that stage, among which Vela X (Blondin et al., 
2001), the Boomerang nebula (Kothes et al., 2006b) and the 
Snail in G327.1-1.1 (Temim et al., 2009, 2015, shown in 
Fig. 10). 

Independently of the pulsar speed, middle aged systems 
are expected to show large asymmetries. The interaction with 
the SNR reverse shock is not expected to happen spherically, 
as simplified one-zone models are forced to assume (see 
Sec. 3.1). The compactness of the PWN contact discontinuity 
itself is partially destroyed by R-T like instabilities (Sec. 2.5) 
well before the onset of reverberation. Of course, in case of a 
high proper motion of the star, the asymmetry is even larger. 
The asymmetry introduced in the PWN morphology in this 
stage is then expected to be a common feature of almost all 
middle aged systems, as well as of old ones, if they have 
not become bow shock nebulae. The original PWN bubble 
might be fragmented, with radio (and gamma-ray) separated 
bubbles surrounding the remaining nebula and expanding 
under adiabatic forces. The background of the PWN can be 
then very noisy, making it difficult to be identified at radio 
or gamma-rays. On the other hand, X-ray emission might 
simply be too faint. 

Estimating the level of fragmentation and mixing of the 
PWN after reverberation can be quite complex: HD simula- 
tions typically do not converge, because in the HD regime 
the fastest growing scale of R-T like instabilities is set purely 
by numerical viscosity at the grid scale (Blondin et al., 2001; 
Bucciantini et al., 2005), while MHD simulations have not 
been performed for old systems; observationally the mix- 
ing has only been estimated for the Snail (Ma et al., 2016), 
suggesting a pulsar wind filling factor of order of 5096. For 
a discussion of the possible impact of mixing on thin-shell 
modelling see Bucciantini et al. (2011). 


5.3 Bow shock nebulae 


Through the combination of radio, H, (in case of a partially 
ionized ambient medium) and, especially, X-ray observations, 
nowadays we have identified 25 fast moving pulsars with an 
associated bow shock nebula. They are listed in Table 2 and 
are marked with blue colored circles in Fig. 8. 

Few bow shock nebulae show an elongated X-ray tail, in 
some cases associated with an even more extended radio tail. 
Only a part of them shows a spectral variation along the X-ray 
tail, in particular a softening indicating synchrotron cooling 
(e.g. the Mouse and the Lighthouse, Kargaltsev et al. 2017). 
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Figure 9. Left panel: Composite optical image of the Crab nebula, Credits: ESO. Right panel: Combined optical (in red — from Hubble) and X-ray (in blue — 
from Chandra) images of the Crab nebula, Credits: Optical - NASA/HST/ASUJ/J. Hester et al. ; X-Ray - NASA/CXC/ASUJJ. Hester et al. 


Figure 10. Composite IR (for the stellar field), radio (in red color) and 
X-ray (in blue) image of the PWN G327.1-1.1, one of the very few sys- 
tems in clear interaction with the SNR reverse shock. Credits: X-Ray — 
NASA/CXC/SAO/T. Temim et al. and ESA/XMM-Newton; Radio: SIFA/MOST 
and CSIRO/ATNF/ATCA; IR: UMass/IPAC-Caltech/NASA/NSF/2MASS. 


To date no TeV emission has been detected from bow shocks 
directly, while UV has only been detected in correspondence 
with Hy emission, probably coming from the heated shocked 
ISM. The bow shock head appears not to have a standard 
morphology, with even drastic variations from one object 
to another. As shown by 3D MHD models (Barkov et al., 
2019; Olmi & Bucciantini, 2019a,b), these variations can be 
ascribed to intrinsic differences in the geometry of the pulsar 
magnetosphere, in the orientation of the pulsar spin-axis with 
respect to the pulsar direction of motion, and the orientation 
with respect to the observer's line of sight. 

In some cases the structure of the bow shock appears to be 
modified in the so called head-and-shoulder shape: the bow 
shock shows an evident widening with distance from the pul- 
sar, with possibly a periodic structure, as the famous example 
of the Guitar nebula (Chatterjee & Cordes, 2004; van Kerk- 
wijk & Ingle, 2008, see e.g. Fig. 11). This is believed to be 
the sign of the mass loading of ambient neutral atoms into the 
bow shock through the shocked ISM (Morlino et al., 2015); 
those atoms then interact with the pulsar wind and modify its 
dynamics (Bucciantini & Bandiera, 2001; Bucciantini, 2002). 
This effect has been proved through numerical simulations by 
Olmi et al. (2018), showing that the lateral expansion of the 
bow shock tail is a function of the pulsar Mach number only, 
namely it increases with the Mach number as the effect of 
the augmented ram pressure exerted by the ISM on the bow 
shock nebula contact discontinuity. 

In recent years bow shock nebulae have gained renewed 
interest thanks to the detection of collimated, extended and 
generally highly misaligned (with respect to the direction of 
motion), jet-like features, only visible at X-rays (observed 
in the Chandra band: 0.5-8 keV), usually referred to as mis- 
aligned tails (De Luca et al., 2011; Pavan et al., 2014; Klin- 
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Figure 11. Composite optical (in blue color) and X-ray (in magenta) image 
of the Guitar BSPWN, generated by the fast moving pulsar J2225+6535, 
and its extended X-ray misaligned tail. The Guitar nebula is also one of 
the systems characterized by the modification of the tail structure in the so 
called head-and-shoulder shape, possibly indicating mass loading from the 
ambient medium into the tail. Credits: X-Ray — NASA/CXC/UMass/S. Johnson 
et al; Optical: NASA/STScI & Palomar Observatory 5-m Hale Telescope. 


gler et al., 2016; de Vries & Romani, 2020; Wang, 2021; de 
Vries et al., 2022; de Vries & Romani, 2022). These structures 
were already observed few years ago surrounding a couple 
of systems (e.g. the Lighthouse nebula), but with the recently 
increased number of detection they now seem a rather com- 
mon feature of these evolved nebulae. At present a generally 
accepted interpretation (Bandiera, 2008) is that they are pro- 
duced by high energy particles (close to the maximum limit 
of the potential drop) leaking from the bow shock nebula, 
and then producing emission via synchrotron radiation in the 
local magnetic field. The observed asymmetry appears to be 
related to the mutual inclination of the spin-axis of the pulsar, 
its magnetic field, the pulsar speed and the direction of the 
magnetic field lines in the ambient medium (Olmi & Buc- 
ciantini, 2019c). Indeed an open question is how to amplify 
the magnetic field from the ISM value to the order of few 
tens of uG required to produce the observed emission. 

Evolved pulsars are also associated with the formation of 
TeV halos (Abeysekara et al., 2017), that have been again 
interpreted as high-energy particles escaping from the bow 
shock nebula, and then diffusing (with some suppression) 
in the ambient medium to form bright gamma-ray bubbles, 
much more extended than the original bow shock (see e.g. a 
sketch for the case of Geminga in Fig. 12). 


6 WHERE WE ARE AND WHERE WE GO 
6.1 Observational prospects 


The actual population of PWNe is expected to increase largely 
in the next future, especially thanks to very high energy (> 
100 GeV) observations. The actual Galactic plane survey 
from the H.E.S.S. telescope (Abdalla et al., 20182) found 
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~30 pc @ 100 TeV 


Figure 12. Sketch roughly comparing the size of the TeV halo around the 
Geminga pulsar wind nebula (from HAWC measures at 100 TeV) and that 
of the X-ray pulsar wind nebula (image adapted from the original composite 
picture at X-rays — Chandra — and IR — Spitzer). Credits for the PWN map: 
X-Ray — NASA/CXC/PSU/B. Posselt et al; IR: NASA/JPL-Caltech. 


that more than half of the 24 extended sources detected are 
identifiable with PWNe (14, mainly thanks to their multi- 
wavelength counterpart). In the Fermi-LAT 3FGL catalog 
(Abdo et al., 2013b) the unidentified sources are ~ 20% of 
the total and it is plausible that most (all?) of them are actually 
PWNe with no direct association with a known pulsar. As 
we already said, thanks to their longer lifetime as gamma- 
ray emitters, middle aged PWNe will likely dominate the 
very high energy sky, possibly representing up to ~ 60% of 
the Galactic sources, and a huge number of new detection 
(^ 200) may be expected with the upcoming Cherenkov 
Telescope Array, thanks to its unprecedented angular and 
energy resolution (Remy et al., 2022; Fiori et al., 2022). A 
very challenging problem will be that of the source confusion 
in the crowded Galactic plane, that will reduce the number of 
identified sources with respect to theoretical estimate (Mestre 
et al., 2022). Moreover an additional source of confusion 
might be that of TeV halos associated with evolved pulsars. 
Given their extended and weak emission, they are difficult to 
identify, and likely constitute a background noise for other 
sources. However the number of expected halos in the Galaxy 
is still matter of debate, ranging from few hundreds to a few, 
depending on their physical interpretation (Sudoh et al., 2019; 
Giacinti et al., 2020; Martin et al., 2022). 

A very exciting recent result is the detection of PeV pho- 
tons coming from 12 sources in the Galaxy, plus the Crab 
nebula, by Cao et al. (2021a,b). Unfortunately, the limited 
angular resolution of LHAASO does not permit to identify 
the exact location of the source (and origin of these energetic 
photons), except for the case of the Crab. One or more pulsars 
can be found in the same region covered by the PSF of the 
instrument for all the 12 sources. Out of these, 11 result to be 
theoretically compatible with being powered by a pulsar (de 
Ofia Wilhelmi et al., 2022), meaning that pulsars and their 
nebulae might also be the most numerous class of extremely 
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high energy emitters in the Galaxy. 

PeV data are also fundamental to either confirm or exclude 
the presence of an hadronic component in the pulsar wind 
(Atoyan & Aharonian, 1996; Bednarek & Protheroe, 1997; 
Bednarek & Bartosik, 2003; Amato et al., 2003). If present, 
hadrons might show up only at the very high energies, where 
the leptonic emission from IC drastically falls due to Klein- 
Nishina suppression (Amato & Olmi, 2021), thus the next 
generation of imaging atmospheric Cherenkov telescopes will 
have a crucial role in answering this question. 

The recently launched IXPE satellite (Weisskopf et al., 
2022), operating in the range 2-8 keV, would enable us to 
sample and image, for the first time, the magnetic field struc- 
ture (not just its strength) and the level of turbulence in a 
handful of PWNe, resolving the magnetic field pattern in 
the central region of the torus-arcs characteristic of young 
objects, and possible in the jets. The three main targets of the 
mission for space resolved polarimetry are: the Crab nebula, 
Vela and MSH 15-52. Other PWNe will be observed in the 
second and third year of operations. Prior to IXPE the Crab 
nebula was the only object to have been observed with a de- 
tected average X-ray polarization of 19%. IXPE will open a 
new observational window into the way we understand and 
characterize these objects. Preliminary results are coming in 
these same days (e.g Bucciantini et al., 2022). By the end of 
the mission we expect to have a much better understanding 
of the dynamical conditions in these relativistic accelerators. 

With the approaching end of operations of the Chandra 
telescope, a very important instrument for future observations 
of PWNe will be the LYNX X-ray observatory$, that thanks 
to its improved sensitivity and field of view promises to open 
new windows of opportunity to investigate the structure and 
details of X-ray sources (an example is the possibility to 
finally detect the compact object in SNR 1987A, see e.g. 
Greco et al. 2021). 


6.2 Modelling prospects 


Despite the astonishing progresses made in the last two 
decades, there are still many important open questions in our 
understanding of the physical processes operating in PWNe. 
If a part of the historical open problems seem nowadays to be 
solved, as the case of the long standing sigma-paradox, some 
are still not properly answered, among which: 


e What are the physical mechanisms responsible for particle 
acceleration at the different energies? AII the proposed 
mechanisms have strengths and weaknesses and require 
precise — and very diverse — physical properties of the 
nebular plasma to be viable. 

What is the origin of the gamma-ray flares in the Crab 
nebula? And is the Crab the unique source producing 
flares? Many possibilities had been proposed and investi- 
gated, most of them requiring a mG magnetic field at the 
emitting region (this is the case for powerful reconnection 
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events in c > 1 regions), much larger than that esti- 
mated from PeV emission. A very detailed discussion of 
this point, and of the proposed models, has been recently 
reported in Amato & Olmi (2021). 

e Are hadrons present in the pulsar wind? Despite being 
certainly a minority by number, hadrons could even be — 
if present — energetically dominant in the wind, changing 
completely our understanding of its properties. 


In the latest years new questions have been added, especially 
thanks to observations of evolved systems that have revealed 
a number of unexpected features: 


e How and with what efficiency particles can escape from 
evolved pulsar wind nebulae? 

How escaped particles produce misaligned tails or ex- 
tended halos? 

Do pulsars accelerate particles at the theoretical limit of 
the maximum potential drop? 

How does the interaction with the reverse shock of the 
SNR modify the emitting properties and morphology of 
evolved PWNe? 

Which ingredients of current modelling need to be mod- 
ified/improved to interpret the upcoming large amount 
of gamma-ray observations and to manage source confu- 
sion? Is there a way to theoretically model the possible 
different morphologies and spectra of evolved PWNe that 
can help in their identification, especially in the lack of a 
multi-wavelength counterpart? 


Recent modelling efforts go in the direction of trying to an- 
swer these questions. How particles escape from evolved 
systems seems now to be clarified, as discussed in Olmi & 
Bucciantini 2019c. The efficiency of the escape process is 
strictly linked to the particle energy, with only the most en- 
ergetic particles, close to the maximum energy achievable, 
able to escape in large fractions. The possibility for them to 
be revealed as extended and diffuse, or asymmetric and thin, 
structures is also partially accounted for with numerical mod- 
els, that predict a variety of escaping processes depending 
on the properties of the system and those of the surrounding 
medium. 

We still lack an understanding of is what happens once 
particles have been injected in the ambient medium. What 
process causes the reduction of the Galactic diffusion length 
in the vicinity of pulsars (if this is the case), to produce 
TeV halos with the observed size through diffusion? And 
also, what amplifies the magnetic field to the value needed 
for the observed synchrotron emission from misaligned X- 
ray tails, a factor of 2 to 10 larger than that expected in the 
ISM? To answer these questions new dedicated modelling 
of particle propagation must be investigated, considering 
diffusion properties, development of self-turbulence and the 
onset of instabilities able to modify the properties of the ISM, 
possibly using a hybrid approach between pure MHD and 
Particle In Cell (PIC) techniques, that correctly accounts for 
the evolution of particles in the plasma. 
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A better understanding of the composition, and of the phys- 
ical properties, of the pulsar wind passes through a correct 
interpretation of high energy gamma-ray data. In particular, 
the firm identification of the observed PeVatrons will shed 
light on the possibility for pulsars to efficiently accelerate par- 
ticles very close to the theoretical limit. A refined modelling 
of the spectral properties through evolutionary phases will 
also be extremely important. 

A very challenging point will be to develop models that 
help in disentangling the PWNe contribution to the gamma- 
ray Galactic emission from that of other sources. Part of this 
requires to understand how to model evolved systems, both in 
terms of emission and morphology, with fast and light enough 
approaches to reproduce a large sample of the expected pop- 
ulation, but enough refined to be reliable. In Bandiera et al. 
(2023) a first step to better include the reverberation phase in 
the description of the PWNe evolution has been made, based 
on the modification of the standard one-zone description. 
Nevertheless these results, despite being much more accurate 
than previous models, still are far for being definitive. A 3D 
MHD modelling of reverberation is necessary to understand 
which role the third spatial dimension, and the structure of the 
magnetic field, play in shaping the PWN during its interaction 
with the SNR. 

Moreover one-zone models by construction cannot account 
for the formation of asymmetric systems, that we expect will 
constitute the largest part of middle aged to evolved PWNe. 
Then, models must be somehow generalized to account for 
different geometries, but how without running expensive 3D 
— or even 2D — models is absolutely unclear, with only few 
preliminary studies presented for the moment (e.g. Olmi & 
Torres 2020). This will also be the flip of the coin for the 
interpretation of gamma-ray data. 


7 CONCLUSIONS 


Pulsar wind nebulae are extremely fascinating systems, show- 
ing a large variety of intriguing properties that require com- 
plex physics to be interpreted. They are known to be powerful 
and efficient particle accelerators and antimatter factories in 
the Galaxy, maybe the primary source of the positron excess 
in the cosmic ray spectrum. They are also the largest class of 
gamma-ray emitting sources in the Galaxy, possibly both at 
very high energies (2 100 GeV) and extremely high energies 
(2 100 TeV). Evolved PWNe are now known to be associated 
with the efficient leakage of particles in the ambient medium, 
showing up in two distinct ways: elongated, thin and asym- 
metric X-ray misaligned tails, originating from the head of 
bow shock PWNe; diffuse and very extended TeV halos, for 
the moment detected around few evolved pulsars. 
Understanding and modelling the properties of PWNe in 
their late evolutionary phases passes through the correct mod- 
elling of all their previous stages. Here we have reviewed 
what we have learned about these different phases, both from 
the observational and theoretical points of view, with focus 
on state of the art numerical modelling. In Section 2 we have 
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in particular posed the bases for the description of the dif- 
ferent stages of a PWN evolution, discussing from a more 
qualitative point of view what characterizes the four phases 
in which we can roughly divide it: the free-expansion phase 
(Sec. 2.3), the reverberation phase (Sec. 2.4), the transitional 
phase after reverberation (Sec. 2.5) and the final phase outside 
the SNR (Sec. 2.6). In Section 2.7 we also reviewed what we 
know about PWNe in other environments, showing how these 
systems are prototypical of many other high-energy astro- 
physical sources. A more quantitative discussion about how 
PWNe through their phases have been modeled can be found 
in Section 3. Here we reviewed all the different approaches 
used up to present days, and highlighted the main results from 
recent 3D MHD numerical simulations of young and old (bow 
shocks) PWNe (see Sec.3.5 and 3.6). A description of the radi- 
ation mechanisms producing the observed emission, and what 
we have understood about the underlying acceleration mecha- 
nisms can be found in Sec. 4. PWNe observational properties, 
and the differences between the various phases, was discussed 
in Section 5. Finally, in Section 6 we present our view about 
observational and theoretical/numerical prospects. 

New impetus was already impressed in last years by the 
observation of TeV halos, misaligned X-ray tails and, very re- 
cently, by the detection of numerous PeVatrons in the Galaxy 
by LHAASO. We believe that many other new clues will 
come with the next generation of IACTs, as CTA or the AS- 
TRI Mini-Array. A very important challenge for the high 
energy astrophysics community will then be the interpreta- 
tion of new gamma-ray data in the coming future: this might 
finally help solving many of the questions that remains unan- 
swered in the fascinating pulsar wind nebulae zoo. 
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A A CATALOG OF PWNE 


In this Appendix we collect the information about all the 
identified Galactic PWNe that we were able to find in the 
literature, considering different observational bands. Here we 
only report those systems for which the association with a 
PWN seems clear, not including the large number of sources 
marked as possible PWNe. Instead a catalog of these systems 
can be found for example in Kargaltsev et al. (2013). 

The source of the various information is specified in the 
caption of the two Tables. In particular we decided to report 
separately PWNe clearly associated with fast moving pulsars 
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(in Table 2) and all the others (in Table 1). 
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Table 1 List of detected PWNe with associated PSR. Values for P, P, B, É and the distance d are taken from the ATNF catalogue, version 1.67. For the X-ray luminosity we report, when 
available, the measure in the 2.1 — 10 keV band (from the Chandra catalog of Galactic sources: hea- www. harvard.edu/ChandraSNR/snrcat. gal.html), otherwise the 0.5 — 8 keV data 
from Kargaltsev et al. (2013, 2017). In some case the luminosity in the 2.1 — 10 keV band is not limited to the PWN and there is a possible contamination from the SNR (marked with a t 
apex, standing for total). If the PWN has been observed in other bands, the information is given in the last column, with: R for radio, O for optical, y for gamma-rays, and H, (data from “the 
Pulsar Wind Nebula Catalog" — www.physics.mcgill.ca/-pulsar/pwncat.html and the TeVCat catalog — tevcat2.uchicago.edu, to which we refer for updated references on the 
instruments detecting the various sources at gamma-rays). A question mark at the apex indicates a non clear association and detection or an uncertain measure. 


SNR PWN PSR P P B É d Lx band also seen 

# S 10-4^ss! 102G ergs! kpc erg sl keV @ 

1  G184.6-05.8 Crab J0534+2200 0.03339 42.10 3.79  4.5x10°8 2 601x109 21-10 R,O,y 
2 G130.7+3.0 3C 58 J0205+6449 0.06572 19.38 3.61  2.7x10?7 3.2 1.10x10%4 — 2.1—10 Ry 
3 G180.0-1.7 G179.72-1.69 J0538+2817 0.1432 0.3669 0.733 49x10? 1.3 1.006x10! 0.5-8 R’ 
4  G2304-01.4! G230.39-1.42 J0729-1448 0.2517 11.33 5.4 2.8x10? 2.7 1.58x10?! 0.5-8 y 
5 G263.9-3.3 Vela J0835-4510 0.08933 12.50 3.38 | 6.9x10?6 0.28 420x102  À2.1—10 R, y 
6 | G284.0-01.8 | G284.08-1.88 71016-5857 0.1074 8.083 2.98 | 2.6x10?6 3.16 1.84x102  2.1— 10° y! 
7 | G2874400.5 Puppy J1048-5832 0.1237 9.612 3.49 2x1036 2.9 844x109! — 2.1 — 10° y 
8  G2922-005  G292.15-0.54 71119-6127 0.408 402.0 41 2.3x1036 8.4 2711x103 2.1 — 10° y! 
9 G292.041.8 | G292.04«41.75  J1124-5916 0.1355 75.25 102. 12x10? 5 1.36x10?5  2.1— 10° Ry’ 
10 G304.1-00.2 G304.10-0.24 71301-6305 0.1845 26.67 7.1 1.7x10°6 10.7 145x107 0.5-8 y 
11  G309.9-02.5 | G309.92-2.51 71357-6429 0.1661 36.02 7.83 3.1x10% 3.1 826x10?! 2.1 — 10° R, y 
12. G313.6+00.3 Kookaburra 71420-6048 0.06818 8.3167 241 1x10°7 5.6 1.4x10%3 0.5-8 R, y 
13 G313.3+00.1 Rabbit J1418-6058 0.11057 16.94 438  4.9x1036 2-5 2.56x10°3 2-10 R, y 
14 | G3204-12 Jellyfish J1513-5908 0.1516 152.9 15.4 17x10? 4.4 1.71x10%5 24-10 R, y 
15  G3324-00.4  G332.50-028 71617-5055 0.06936 13.51 3.1 1.6x10?7 4.7 6.08x10?  2.1—10 - 
16 | G344.7-00.1 G344.7440.12 J1702-4128 0.1821 5.234 3.12 3.45105 4.0 4.84x10?5  2.1— 10° y! 
17 G348.9-00.4"  G348.95-0.43 71718-3825 0.07467 1.322 1.00 — L3x10?6 3.5 3.98107? 0.5-8 - 
18 G034.04202 G34.01+20.27 J1740+1000 0.1541 2.147 1.84 2.3x1035 1.2 1.21x102! — 2.1 — 10° R 
19  G008.3400.] | G8.440«40.15 71803-2137 0.1337 13.44 429  22x10*6 44 2.8107 — 2.1 — 10! y 
20 GO011.1+00.1 — G11.0940.08 — 71809-1917 0.08276 2.553 147 1.8x10°6 3.3 458x107 241-10 — Ry 
21  GO011.2-00.3 Turtle J1811-1925 0.06467 4.400 1.731 | 64x10?6 5.0 894x109 2.1 — 10° Ry’ 
22  G018.0-00.6  G18.00-0.69 71826-1334 0.1015 7.525 2.8 | 2.8x10?6 3.6 4.69102  2.1—10 R, y 
23 | G21.5-0.9 G21.50-0.89 . 71833-1034 0.06188 20.20 3.38 | 3.4x10?? 4.1 2.9x105  2.1—10 R, y 
24  G029.7-00.2 Kes75 J1846-0258 0.3266 710.7 48.8 | 8.1x10?6 5.8 140x106  2.1—10 R, y 
25  G034.7-04 G34.56-0.50 7185640113 0.2674 20.84 7.55 | 43x10 3.3 497x102  2.1—10 R 
26  G054.140.3 | G54.10+0.27 71930-1852 0.1369 75.06 103 12x10? 7 2.03x10°* — 2.1— 10° Ry 
27  G047.3-03.8 | G47.38-3.88 7193241059 0.2265 0.1157 0.518 — 3.910? 0.31 9.89x1022 — 2.1— 10° R’ 
28 G075.2+00.1 Dragonfly J2021+3651 0.1037 9.572 3.19 3.4x10?6 1.8 741x102? 2.1-10 R-y’ 
29 G106.6+02.9 Boomerang J2229+6114 0.05162 7.827 2.0 22x10? 3.0 73x102 2. -10 Ry’ 
30  G012.8-0.00 | G12.82-0.02 71813-1749 0.04474 12.70 241 5.6x1037 6.15 5.35x10°* 2.1 — 10° y 
31 G310.6-1.6 G310.6-1.6 J1400-6325 0.03118 3.890 1.11 — 5.1x10? 7.0 11340? 2.1 — 10° R 
32 SNR W42 G25.24-0.19 71838-0655 0.0705 4.925 1.89 — 5.5x10?6 6.6 = E y 
33 G119.5+10.2 CTAI J00074-7303 0.3159 36.00 10.8 | 4.5x10?5 1.5 245x10! — 2.1 — 10° y 
34  G007.5-01.7? Taz J1809-2332 0.1468 3.442 227  43x105 0.88?! 127x109  2.1-— 10° R 
35  G266.9-01.0 | G266.97-1.00  J0855-4644 0.06469 0.7263 0.694 1.1x10° 0.5-5.6? L04x102  2.1— 10° R 
36 G021.9-00.1? | G21.88-0.1 J1831-0952 0.06727 . 0.8324 0.757  lL.bxd0?6 3.7 = = y 
37  G076.9-01.0 = J2022+3842 0.04858 8.610 2.07  3.0x107 — 7-107 442x10? — 2.1 — 10° R 


vc 


nuupioong `N 3 nuo `J 


Table 2 List of the known pulsars with high proper motion and an associated PWN. Pulsar values are taken, as before, from ATNF catalogue, version 1.67, while the X-ray luminosity is 
always taken from Kargaltsev et al. (2013) and Kargaltsev et al. (2017). Symbols and notation are the same defined in the previous table. 


PWN/ASSOCIATED OBJ PSR P P B É d Lx band also seen 
# s 10714 ss! 10! G ergs! kpc erg s! keV @ 
1 Geminga J0633+1746 0.2371 1.097 1.63 32x10^ 0.19 224x10? 0.5-8 y 
2 G319.97-0.62 71509-5850 . 0.08892 0.9166 0.914 5.x10? 34  L12x10? 0.5-8 R 
3 G343.10-2.69 J1709-4429 — 0.1025 9.298 3.12 3.4x106 2.6 398xl10? 0.5-8 Ry 
4 Mouse J1747-2958 0.09881 6.132 2.49 2.5106 2.5 6.76x10 05-8 R 
5 Duck J1801-2451 0.1249 12.79 4.04 2.6x10°° 3.8  L6xl0?5 05-8 R, y 
6 Guitar J222546535 0.6825 0.9661 2.6 1.22340? 0.8%  L51x10?9 0.5-8 S 
7 Morla J035743205 0.4441 1.304 2.43 59x10? 0.8  L17x10?9? 0.5-8 = 
8 - J0437-4715 0.005757  5.7x105  5.81x10^ 1.2x10^ 0.16 ~~ 4x108 05-8 Ha, R? 
9 SNR S147 J053842817 . 0.1432 0.3669 0.733 49xl0^ 1.3 2x10?! 0.5-8 
10 E J0742-2822 0.1668 1.682 1.69 1.4x105 2 - = H,, R? 
11 - J0908-4913 . 0.1068 1.51 1.28 49xl0^ 1 - - H?,R 
12 Lighthouse J1101-6101 0.0628 0.86 0.742 1.4x10°° 7 2.x10? 0.5-8 H? 
13 G293.79+0.58 J1135-6055 0.1149 7.93 3.05 2.1x10% 29 25x10? 05-8 HLR 
14 Frying Pan J1437-5959 0.0617 0.8587 0.737 1.4x1076 — 8,5 - = R 
15 G006.4+04.9 J1741-2054 0.4137 1.698 2.68 9.510? 03  L6xl0? 05-8 2 
16 Eel J1826-1256 0.1102 12.15 3.7 3.6106 16 2.4x10% 05-8 HLy 
17 SNR W44 J1856+0113 — 0.2674 20.84 7.55 43xl0?5 3.3 16x10" 05-8 R 
18 G47.38-3.88 J1932+1059 — 0.2265 0.1157 0.518 39x10? 03  Á3.16x10? 0.5-8  HLR 
19 SNR CTB 80 J1952+3252 0.03953 0.5845 0.486 3.7x10°° 3 1x10? 05-8 HR 
20 Black Widow J1959+2048 0.001607  1.69x10-.6  1.67x10-^ 1.6x10%5 14  5.5d0? 05-8 R 
21 - J2030+4415 02271 0.6484 1.23 22x10^ 0.7  3.x10? 05-8 H, 
22 - J205542539 0.3196 0.408 1.16 4.9x10% 06 1.48x10 05-8 HLR' 
23 G10.92-45.43 J2124-3358 0.004931 2.06x10-6 322x107 6.8x10? 0.4  9.5x108 05-8 R 
24 Mushroom J0358+5413 . 0.1564 0.4395 0.839 4.5x10°4 1 1.58x10! 0.5-8 - 
25 = J1648-4611 0.165 2.373 2.0 2.00? 45 <3xl0! 05-8 = 


^ Updated distance from Deller et al. (2019). 


INMd Jo ymd &ipnuonnjoaa ay J 


SC 


26 


REFERENCES 


Abdalla H., Abramowski A., Aharonian F., et al. for the 
H. E. S. S. Collaboration 2018a, A&A, 612, A1 

Abdalla H., Abramowski A., Aharonian F., et al. for the 
H. E. S. S. Collaboration 2018b, A&A, 612, A2 

Abdo A. A., et al., 2013a, ApJS, 208, 17 

Abdo A. A., et al., 2013b, ApJS, 208, 17 

Abeysekara A. U., et al., 2017, Science, 358, 911 

Actis M., Agnetta G., Aharonian e. a., 2011, Experimental 
Astronomy, 32, 193 

Aharonian F., et al., 2005, A&A, 435, L17 

Amato E., 2014, in International Journal of Modern Physics 
Conference Series. p. 1460160 (arXiv:1312.5945), 
doi:10.1142/82010194514601604 

Amato E., Olmi B., 2021, Universe, 7, 448 

Amato E., Guetta D., Blasi P., 2003, A&A, 402, 827 

Archibald R. F., et al., 2016, ApJ, 819, L16 

Arons J., 2012, Space Sci. Rev., 173, 341 

Arzoumanian Z., Chernoff D. F., Cordes J. M., 2002, ApJ, 
568, 289 

Atoyan A. M., 1999, A&A, 346, L49 

Atoyan A. M., Aharonian F. A., 1996, MNRAS, 278, 525 

Aumont J., et al., 2010, A&A, 514, A70 

Baade W., 1956, Bull. Astron. Inst. Netherlands, 12, 312 

Bandiera R., 2008, A&A, 490, L3 

Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 
2020, MNRAS, 499, 2051 

Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 
2021, MNRAS, 508, 3194 

Bandiera R., Bucciantini N., Martín J., Olmi B., Torres D. F., 
2023, MNRAS, 

Barkov M. V., Komissarov S. S., 2008, International Journal 
of Modern Physics D, 17, 1669 

Barkov M. V., Lyutikov M., Khangulyan D., 2019, MNRAS, 
484, 4760 

Bednarek W., Bartosik M., 2003, A&A, 405, 689 

Bednarek W., Protheroe R. J., 1997, Phys. Rev. Lett., 79, 
2616 

Bednarek W., Sitarek J., 2013, MNRAS, 430, 2951 

Begelman M. C., 1998, ApJ, 493, 291 

Begelman M. C., Li Z.-Y., 1992, ApJ, 397, 187 

Bietenholz M. F., Kronberg P. P., 1990, ApJ, 357, L13 

Bietenholz M. F., Frail D. A., Hankins T. H., 1991, ApJ, 376, 
L41 

Bietenholz M. F., Frail D. A., Hester J. J., 2001, ApJ, 560, 
254 

Bietenholz M. F., Hester J. J., Frail D. A., Bartel N., 2004, 
ApJ, 615, 794 

Blondin J. M., Chevalier R. A., Frierson D. M., 2001, ApJ, 
563, 806 

Blumer H., Safi-Harb S., McLaughlin M. A., 2017, ApJ, 850, 
L18 

Bogovalov S. V., 1999, A&A, 349, 1017 


B. Olmi & N. Bucciantini 


Bogovalov S. V., Khangoulian D. V., 2002, MNRAS, 336, 
L53 

Bogovalov S. V., Khangoulyan D. V., 2002, Astronomy Let- 
ters, 28, 373 

Bogovalov S. V., Chechetkin V. M., Koldoba A. V., Ustyu- 
gova G. V., 2005, MNRAS, 358, 705 

Bosch-Ramon V., Barkov M. V., Khangulyan D., Perucho M., 
2012, A&A, 544, A59 

Bucciantini N., 2002, A&A, 387, 1066 

Bucciantini N., 2018, MNRAS, 478, 2074 

Bucciantini N., Bandiera R., 2001, A&A, 375, 1032 

Bucciantini N., Blondin J. M., Del Zanna L., Amato E., 2003, 
A&A, 405, 617 

Bucciantini N., del Zanna L., Amato E., Volpi D., 2005, A&A, 
443, 519 

Bucciantini N., Quataert E., Arons J., Metzger B. D., Thomp- 
son T. A., 2007, MNRAS, 380, 1541 

Bucciantini N., Quataert E., Arons J., Metzger B. D., Thomp- 
son T. A., 2008, MNRAS, 383, L25 

Bucciantini N., Arons J., Amato E., 2011, MNRAS, 410, 381 

Bucciantini N., Metzger B. D., Thompson T. A., Quataert E., 
2012, MNRAS, 419, 1537 

Bucciantini N., et al., 
arXiv:2207.05573 

Camero-Arranz A., et al., 2013, MNRAS, 429, 2493 

Camilo F., Gaensler B. M., Gotthelf E. V., Halpern J. P., 
Manchester R. N., 2004, ApJ, 616, 1118 

Camus N. F., Komissarov S. S., Bucciantini N., Hughes P. A., 
2009, MNRAS, 400, 1241 

Cao Z., et al., 2021a, Science, 373, 425 

Cao Z., et al., 2021b, Nature, 594, 33 

Cerutti B., Philippov A. A., 2017, A&A, 607, A134 

Cerutti B., Philippov A. A., Dubus G., 2020, A&A, 642, 
A204 

Chatterjee S., Cordes J. M., 2004, ApJ, 600, L51 

Chen L., Zhang B., 2021, ApJ, 906, 105 

Chevalier R. A., 1976, ApJ, 207, 872 

Childs H., et al., 2012 

Ciolfi R., Siegel D. M., 2015, ApJ, 798, L36 

Contopoulos I., Kazanas D., Fendt C., 1999, ApJ, 511, 351 

Cordes J. M., Chernoff D. F., 1998, ApJ, 505, 315 

De Luca A., et al., 2011, ApJ, 733, 104 

Del Zanna L., Amato E., Bucciantini N., 2004, A&A, 421, 
1063 

Del Zanna L., Volpi D., Amato E., Bucciantini N., 2006, 
A&A, 453, 621 

Deller A. T., et al., 2019, ApJ, 875, 100 

Dessart L., 2019, A&A, 621, A141 

Dessart L., Hillier D. J., Waldman R., Livne E., Blondin S., 
2012, MNRAS, 426, L76 

Dodson R., Lewis D., McConnell D., Deshpande A. A., 2003, 
MNRAS, 343, 116 

Emmering R. T., Chevalier R. A., 1987, ApJ, 321, 334 


2022, arXiv e-prints, p. 


The evolutionary path of PWNe 


Faucher-Giguére C.-A., Kaspi V. M., 2006, ApJ, 643, 332 

Ferreira S. E. S., de Jager O. C., 2008, A&A, 478, 17 

Fiori M., Zampieri L., Burtovoi A., Caraveo P., Tibaldo L., 
2020, MNRAS, 499, 3494 

Fiori M., Olmi B., Amato E., Bandiera R., Bucciantini N., 
Zampieri L., Burtovoi A., 2022, MNRAS, 511, 1439 

Frail D. A., Giacani E. B., Goss W. M., Dubner G., 1996, 
ApJ, 464, L165 

Gaensler B. M., Slane P. O., 2006, ARA&A, 44, 17 

Gaensler B. M., Brazier K. T. S., Manchester R. N., Johnston 
S., Green A. J., 1999, MNRAS, 305, 724 

Gaensler B. M., Arons J., Kaspi V. M., Pivovaroff M. J., 
Kawai N., Tamura K., 2002, ApJ, 569, 878 

Gaensler B. M., et al., 2005, Nature, 434, 1104 

Gelfand J. D., 2017, in Torres D. F., ed., Astrophysics and 
Space Science Library Vol. 446, Modelling Pulsar Wind 
Nebulae. p. 161, doi:10.1007/978-3-319-63031-1 8 

Gelfand J. D., et al., 2005, ApJ, 634, L89 

Gelfand J. D., Slane P. O., Zhang W., 2009, ApJ, 703, 2051 

Giacinti G., Mitchell A. M. W., López-Coto R., Joshi V., 
Parsons R. D., Hinton J. A., 2020, A&A, 636, A113 

Goldreich P., Julian W. H., 1969, ApJ, 157, 869 

Gompertz B. P., O'Brien P. T., Wynn G. A., 2014, MNRAS, 
438, 240 

Granot J., Gill R., Younes G., Gelfand J., Harding A., Kouve- 
liotou C., Baring M. G., 2017, MNRAS, 464, 4895 

Greco E., et al., 2021, ApJ, 908, L45 

Guest B., Safi-Harb S., MacMaster A., Kothes R., Olmi B., 
Amato E., Bucciantini N., Arzoumanian Z., 2020, MNRAS, 
491, 3013 

Helfand D. J., Gotthelf E. V., Halpern J. P., 2001, ApJ, 556, 
380 

Hester J. J., 2008, ARA&A, 46, 127 

Hester J. J., et al., 1995, ApJ, 448, 240 

Hester J. J., et al., 2002, ApJ, 577, L49 

Horvath J. E., 2019, MNRAS, 484, 1983 

Huber D., Kissmann R., Reimer A., Reimer O., 2021a, A&A, 
646, A91 

Huber D., Kissmann R., Reimer O., 2021b, A&A, 649, A71 

Jun B.-I., 1998, ApJ, 499, 282 

Kagan D., Sironi L., Cerutti B., Giannios D., 2015, 
Space Sci. Rev., 191, 545 

Kargaltsev O., Pavlov G. G., 2008, in Bassa C., Wang Z., 
Cumming A., Kaspi V. M., eds, American Institute of 
Physics Conference Series Vol. 983, 40 Years of Pulsars: 
Millisecond Pulsars, Magnetars and More. pp 171—185 
(arXiv: 0801.2602), doi:10.1063/1.2900138 

Kargaltsev O., Rangelov B., Pavlov G., 2013, in , The 
Universe Evolution: Astrophysical and Nuclear Aspects. 
Edited by I. Strakovsky and L. Blokhintsev. Nova Science 
Publishers. pp 359—406 

Kargaltsev O., Pavlov G. G., Durant M., Volkov I., Hare J., 
2014, ApJ, 784, 124 

Kargaltsev O., Pavlov G. G., Klingler N., Rangelov B., 2017, 


27 


Journal of Plasma Physics, 83, 635830501 

Kennel C. F., Coroniti F. V., 1984a, ApJ, 283, 694 

Kennel C. F., Coroniti F. V., 1984b, ApJ, 283, 710 

Khangoulian D. V., Bogovalov S. V., 2003, Astronomy Let- 
ters, 29, 495 

Khangulyan D., Arakawa M., Aharonian F., 2020, MNRAS, 
49], 3217 

Kirk J. G., 2006, Advances in Space Research, 37, 1970 

Kirk J. G., Skjzeraasen O., 2003, ApJ, 591, 366 

Kirk J. G., Lyubarsky Y., Petri J., 2009, The Theory of Pulsar 
Winds and Nebulae. Springer Berlin Heidelberg, Berlin, 
Heidelberg, pp 421—450, doi:10.1007/978-3-540-76965- 
1 16 

Klepser S., et al., 2013, in International Cosmic Ray Confer- 
ence. p. 971 (arXiv:1307.7985) 

Klingler N., et al., 2016, ApJ, 833, 253 

Kolb C., Blondin J., Slane P., Temim T., 2017, ApJ, 844, 1 

Komissarov S. S., 2006, MNRAS, 367, 19 

Komissarov S. S., Lyubarsky Y. E., 2004, MNRAS, 349, 779 

Kothes R., Reich W., Uyaniker B., 2006a, ApJ, 638, 225 

Kothes R., Reich W., Uyaniker B., 2006b, ApJ, 638, 225 

Kothes R., Landecker T. L., Reich W., Safi-Harb S., Arzou- 
manian Z., 2008, ApJ, 687, 516 

Kurono Y., Morita K.-I., Kamazaki T., 2009, Publications of 
the Astronomical Society of Japan, 61, 873 

Lai P. C. W., Ng C. Y., Bucciantini N., 2022, ApJ, 930, 1 

Legg M. P. C., Westfold K. C., 1968, ApJ, 154, 499 

Leung W. Y., Ng C. Y., 2016, in Supernova Remnants: An 
Odyssey in Space after Stellar Death. p. 53 

Li Z.-Y., 1993, PhD thesis, University of Colorado, Boulder 

Lou Y.-Q., 1998, MNRAS, 294, 443 

Lu F. J., Wang Q. D., Aschenbach B., Durouchoux P., Song 
L. M., 2002, ApJ, 568, L49 

Lundqvist N., Lundqvist P., Bjórnsson C. I., Olofsson G., 
Pires S., Shibanov Y. A., Zyuzin D. A., 2011, MNRAS, 
413, 611 

Lyne A. G., Jordan C. A., Graham-Smith F., Espinoza C. M., 
Stappers B. W., Weltevrede P., 2015, MNRAS, 446, 857 

Lyubarskij Y. E., 1992, Soviet Astronomy Letters, 18, 356 

Lyubarsky Y. E., 2002, MNRAS, 329, L34 

Lyubarsky Y. E., 2003, MNRAS, 345, 153 

Lyubarsky Y., 2010, ApJ, 725, L234 

Lyutikov M., Komissarov S. S., Porth O., 2016, MNRAS, 
456, 286 

Lyutikov M., Temim T., Komissarov S., Slane P., Sironi L., 
Comisso L., 2019, MNRAS, 489, 2403 

Ma Y. K., Ng C. Y., Bucciantini N., Slane P. O., Gaensler 
B. M., Temim T., 2016, ApJ, 820, 100 

Manchester R. N., Hobbs G. B., Teoh A., Hobbs M., 2005, 
AJ, 129, 1993 

Margalit B., Metzger B. D., 2018, ApJ, 868, L4 

Margalit B., Metzger B. D., Thompson T. A., Nicholl M., 
Sukhbold T., 2018, MNRAS, 475, 2659 


28 


Martín J., Torres D. F., Rea N., 2012, MNRAS, 427, 415 

Martin P., Tibaldo L., Marcowith A., Abdollahi S., 2022, 
A&A, 666, A7 

Melatos A., 1998, Mem. Soc. Astron. Italiana, 69, 1009 

Mestre E., Torres D. F., de Oña Wilhelmi E., Martí J., 2022, 
MNRAS, 517, 3550 

Metzger B. D., Piro A. L., 2014, MNRAS, 439, 3916 

Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., 
Quataert E., 2011a, MNRAS, 413, 2031 

Metzger B. D., Giannios D., Thompson T. A., Bucciantini N., 
Quataert E., 2011b, MNRAS, 413, 2031 

Mignone A., Striani E., Tavani M., Ferrari A., 2013, MNRAS, 
436, 1102 

Milisavljevic D., Patnaude D. J., Chevalier R. A., Raymond 
J. C., Fesen R. A., Margutti R., Conner B., Banovetz J., 
2018, ApJ, 864, L36 

Moran P., Shearer A., Mignani R. P., Stowikowska A., De 
Luca A., Gouiffés C., Laurent P., 2013, MNRAS, 433, 
2564 

Mori K., Burrows D. N., Hester J. J., Pavlov G. G., Shibata 
S., Tsunemi H., 2004, ApJ, 609, 186 

Moriya T. J., Chen T.-W., Langer N., 2017, ApJ, 835, 177 

Morlino G., Lyutikov M., Vorster M., 2015, MNRAS, 454, 
3886 

Nakamori T., et al., 2008, ApJ, 677, 297 

Nalewajko K., Begelman M. C., 2012, MNRAS, 427, 2480 

Neronov A., Chernyakova M., 2007, arXiv e-prints, pp astro— 
ph/0701144 

Ng C. Y., Romani R. W., 2004, ApJ, 601, 479 

Ng C. Y., Gaensler B. M., Chatterjee S., Johnston S., 2010, 
ApJ, 712, 596 

Olmi B., Bucciantini N., 2019a, MNRAS, 484, 5755 

Olmi B., Bucciantini N., 2019b, MNRAS, 488, 5690 

Olmi B., Bucciantini N., 2019c, MNRAS, 490, 3608 

Olmi B., Torres D. F., 2020, MNRAS, 494, 4357 

Olmi B., Del Zanna L., Amato E., Bandiera R., Bucciantini 
N., 2014, MNRAS, 438, 1518 

Olmi B., Del Zanna L., Amato E., Bucciantini N., 2015, 
MNRAS, 449, 3149 

Olmi B., Del Zanna L., Amato E., Bucciantini N., Mignone 
A., 2016, Journal of Plasma Physics, 82, 635820601 

Olmi B., Bucciantini N., Morlino G., 2018, MNRAS, 481, 
3394 

Oort J. H., Walraven T., 1956, Bull. Astron. Inst. Netherlands, 
12, 285 

Pacini F., Salvati M., 1973, ApJ, 186, 249 

Paredes-Fortuny X., Bosch-Ramon V., Perucho M., Ribó M., 
2015, A&A, 574, A77 

Parthasarathy A., et al., 2020, MNRAS, 494, 2012 

Pavan L., et al., 2014, A&A, 562, A122 

Pavlov G. G., Kargaltsev O. Y., Sanwal D., Garmire G. P., 
2001, ApJ, 554, L189 

Pavlov G. G., Teter M. A., Kargaltsev O., Sanwal D., 2003, 


B. Olmi & N. Bucciantini 


ApJ, 591, 1157 

Petre R., Kuntz K. D., Shelton R. L., 2002, ApJ, 579, 404 

Philippov A. A., Spitkovsky A., 2014, ApJ, 785, L33 

Porth O., Komissarov S. S., Keppens R., 2013, MNRAS, 431, 
L48 

Porth O., Komissarov S. S., Keppens R., 2014, MNRAS, 438, 
278 

Rees M. J., Gunn J. E., 1974, MNRAS, 167, 1 

Remy Q., Tibaldo L., Acero E, Fiori M., Knódlseder J., 
Olmi B., Sharma P., 2022, in 37th International Cos- 
mic Ray Conference. 12-23 July 2021. Berlin. p. 886 
(arXiv: 2109.03729) 

Reynolds S. P., Chevalier R. A., 1984, ApJ, 278, 630 

Rezzolla L., Kumar P., 2015, ApJ, 802, 95 

Roberts M. S. E., Romani R. W., Johnston S., Green A. J., 
1999, ApJ, 515, 712 

Roberts M. S. E., Romani R. W., Johnston S., 2001, ApJ, 561, 
L187 

Romani R. W., Ng C.-Y., 2003, ApJ, 585, L41 

Romani R. W., Ng C.-Y., Dodson R., Brisken W., 2005, ApJ, 
631, 480 

Rowlinson A., O’Brien P. T., Metzger B. D., Tanvir N. R., 
Levan A. J., 2013, MNRAS, 430, 1061 

Rybicki G. B., Lightman A. P., 1979, Radiative processes in 
astrophysics 

Safi-Harb S., Ferrand G., Matheson H., 2013, in van Leeuwen 
J., ed., Vol. 291, Neutron Stars and Pulsars: Chal- 
lenges and Opportunities after 80 years. pp 483—485 
(arXiv:1210.5264), doi:10.1017/S 1743921312024593 

Sartore N., Ripamonti E., Treves A., Turolla R., 2010, A&A, 
510, A23 

Scargle J. D., 1969, ApJ, 156, 401 

Schweizer T., Bucciantini N., Idec W., Nilsson K., Tennant 
A., Weisskopf M. C., Zanin R., 2013, MNRAS, 433, 3325 

Scuderi S., et al., 2022, Journal of High Energy Astrophysics, 
35, 52 

Shankar S., Mösta P., Barnes J., Duffell P. C., Kasen D., 2021, 
MNRAS, 508, 5390 

Siegel D. M., Ciolfi R., 2016a, ApJ, 819, 14 

Siegel D. M., Ciolfi R., 2016b, ApJ, 819, 15 

Sironi L., Spitkovsky A., 2011, ApJ, 741, 39 

Sironi L., Keshet U., Lemoine M., 2015, Space Sci. Rev., 191, 
519 

Slane P., 2017, in Alsabti A. W., Murdin P., eds, , Handbook 
of Supernovae. p. 2159, doi:10.1007/978-3-319-21846- 
5 95 

Slane P., Helfand D. J., van der Swaluw E., Murray S. S., 
2004, ApJ, 616, 403 

Smartt S. J., 2009, ARA &A, 47, 63 

Soker N., Gilkis A., 2017, ApJ, 851, 95 

Spitkovsky A., 2006, ApJ, 648, L51 

Sudoh T., Linden T., Beacom J. F., 2019, Phys. Rev. D, 100, 
043016 


The evolutionary path of PWNe 


Swartz D. A., et al., 2015, ApJ, 808, 84 

Takamoto M., Inoue T., Inutsuka S.-i., 2012, ApJ, 755, 76 

Tanaka S. J., 2016, ApJ, 827, 135 

Tanaka S. J., Takahara F., 2010, ApJ, 715, 1248 

Tanaka S. J., Takahara F., 2013, MNRAS, 429, 2945 

Temim T., Slane P., Gaensler B. M., Hughes J. P., Van Der 
Swaluw E., 2009, ApJ, 691, 895 

Temim T., Slane P., Kolb C., Blondin J., Hughes J. P., Buc- 
ciantini N., 2015, ApJ, 808, 100 

Thompson C., 1994, MNRAS, 270, 480 

Thompson T. A., 2007, in Revista Mexicana de 
Astronomia y Astrofisica, vol. 27. pp 80-90 
(arXiv:astro-ph/0611368) 

Timokhin A. N., 2006, MNRAS, 368, 1055 

Timokhin A. N., Harding A. K., 2019, ApJ, 871, 12 

Toropina O. D., Romanova M. M., Toropin Y. M., Lovelace 
R. V. E., 2001, ApJ, 561, 964 

Toropina O. D., Romanova M. M., Lovelace R. V. E., 2019, 
MNRAS, 484, 1475 

Torres D. F., 2017, ApJ, 835, 54 

Torres D. E., 2018, in Weltevrede P., Perera B. B. P., Preston 
L. L., Sanidas S., eds, Vol. 337, Pulsar Astrophysics 
the Next Fifty Years. pp 255-258 (arXiv:1710.07026), 
doi:10.1017/S1743921317008213 

Torres D. F., Cillis A., Martin J., de Oña Wilhelmi E., 2014, 
Journal of High Energy Astrophysics, 1,31 

Torres D. F., Lin T., Coti Zelati F., 2019, MNRAS, 486, 1019 

Truelove J. K., McKee C. F., 1999, ApJS, 120, 299 

Usov V. V., 1992, Nature, 357, 472 

Uzdensky D. A., 2016, in Gonzalez W., Parker E., eds, 
Astrophysics and Space Science Library Vol. 427, 
Magnetic Reconnection: Concepts and Applications. 
p. 473 (arXiv:1510.05397), doi:10.1007/978-3-319- 
26432-5 12 

Uzdensky D. A., 2018, MNRAS, 477, 2849 

Velusamy T., 1985, MNRAS, 212, 359 

Verbunt F., Igoshev A., Cator E., 2017, A&A, 608, A57 

Vigelius M., Melatos A., Chatterjee S., Gaensler B. M., 
Ghavamian P., 2007, MNRAS, 374, 793 

Vink J., Bamba A., 2009, ApJ, 707, L148 

Volpi D., Del Zanna L., Amato E., Bucciantini N., 2008, 
A&A, 485, 337 

Vurm I., Metzger B. D., 2021, ApJ, 917, 77 

Wang L. J., 2017, Acta Astronomica Sinica, 58, 29 

Wang Q. D., 2021, Research Notes of the American Astro- 
nomical Society, 5, 5 

Wang S.-Q., Wang L.-J., Dai Z.-G., 2019, Research in As- 
tronomy and Astrophysics, 19, 063 

Weisskopf M. C., Silver E. H., Kestenbaum H. L., Long K. S., 
Novick R., 1978, ApJ, 220, L117 

Weisskopf M. C., et al., 2000, ApJ, 536, L81 

Weisskopf M. C., et al., 2022, Journal of Astronomical Tele- 
scopes, Instruments, and Systems, 8, 026002 


29 


Westfold K. C., 1959, ApJ, 130, 241 

Wilkin F. P., 1996, ApJ, 459, L31 

Woltjer L., 1958, PhD thesis, Leiden Observatory 

Zajczyk A., et al., 2012, A&A, 542, A12 

Zdziarski A. A., Neronov A., Chernyakova M., 2010, MN- 
RAS, 403, 1873 

Zhu B.-T., Zhang L., Fang J., 2019, ApJ, 873, 120 

de Jager O. C., Ferreira S. E. S., Djannati-Atai A., 2008, in 
Aharonian F. A., Hofmann W., Rieger F., eds, American 
Institute of Physics Conference Series Vol. 1085, Amer- 
ican Institute of Physics Conference Series. pp 199—202, 
doi:10.1063/1.3076638 

de Ofia-Wilhelmi E., et al., 2013, Astroparticle Physics, 43, 
287 

de Ofia Wilhelmi E., López-Coto R., Amato E., Aharonian 
F., 2022, ApJ, 930, L2 

de Vries M., Romani R. W., 2020, ApJ, 896, L7 

de Vries M., Romani R. W., 2022, ApJ, 928, 39 

de Vries M., et al., 2022, arXiv e-prints, p. arXiv:2210.01228 

van Kerkwijk M. H., Ingle A., 2008, ApJ, 683, L159 

van der Swaluw E., Achterberg A., Gallant Y. A., Tóth G., 
2001, A&A, 380, 309 

van der Swaluw E., Achterberg A., Gallant Y. A., Downes 
T. P., Keppens R., 2003, A&A, 397, 913 

van der Swaluw E., Downes T. P., Keegan R., 2004, A&A, 
420, 937 


